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ABSTRACT 

We present physical results for a variety of light hadronic quantities obtained via a combined 
analysis of three 2+1 flavour domain wall fermion ensemble sets. For two of our ensemble sets we 
used the Iwasaki gauge action with = 2.13 (a" 1 = 1.75(4) GeV) and j8 = 2.25 (a" 1 = 2.31(4) 
GeV) and lattice sizes of 24 3 x 64 and 32 3 x 64 respectively, with unitary pion masses in the 
range 293(5)-417(10) MeV. The extent L s for the 5 th dimension of the domain wall fermion 
formulation is L s = 16 in these ensembles. In this analysis we include a third ensemble set 
that makes use of the novel Iwasaki+DSDR (Dislocation Suppressing Determinant Ratio) gauge 
action at /3 = 1.75 (a -1 = 1.37(1) GeV) with a lattice size of 32 3 x 64 and L s = 32 to reach 
down to partially-quenched pion masses as low as 143(1) MeV and a unitary pion mass of 
171(1) MeV, while retaining good chiral symmetry and topological tunneling. We demon- 
strate a significant improvement in our control over the chiral extrapolation, resulting in much 
improved continuum predictions for the above quantities. The main results of this analysis in- 
clude the pion and kaon decay constants, f n = 127(3) st at(3)sys MeV and fg = 152(3) s tat(2) sys 
MeV respectively (fx/fx = 1.199(12) s tat(14) sys ); the average up/down quark mass and the 
strange-quark mass in the Ms-scheme at 3 GeV, w^ms, 3 GeV) = 3.05(8) s tat(6) S y S MeV and 
m s (Ms,3 GeV) = 83.5(1.7) sta t(l-l) S y S ; the neutral kaon mixing parameter in the Ms-scheme at 3 
GeV, B k (ms, 3 GeV) = 0.535(8) stat (13) sys , and in the RGI scheme, B K = 0.758(1 l) stat (19) sys ; and 
the Sommer scales r x = 0.323(8) stat (4) sys fm and r Q = 0.480(10) stat (4) sys (n/r = 0.673(1 l) stat (3) sys ). 
We also obtain values for the SU(2) chiral perturbation theory effective couplings, I3 = 2.91 (23) sta t(7) S} 
and/"4 = 3.99(16) sta t(9) sys . 
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I. INTRODUCTION 

The RBC & UKQCD collaborations have recently published continuum limit results yj, |2J] for a 
variety of light hadronic quantities, including the pion and kaon decay constants, quark masses 
and the neutral kaon mixing parameter Bk, determined using two ensemble sets of 2 + 1-flavor 
domain wall fermions (DWF) with the Iwasaki gauge action at /3 = 2.25 (corresponding to a lattice 
spacing of a « 0.086 fm) and /3 = 2. 13 (a « 0. 1 14 fm), with lattice sizes of 32 3 x 64 and 24 3 x 64 
respectively and fifth-dimensional extents of L s = 16. We refer to this as the '2010 analysis'. With 
precise non-perturbative renormalization methods made possible by the good chiral symmetry of 
the action, and a combined chiral/continuum fit analysis to maximise the use of the available data, 
our predictions were limited mainly by the €?(5%) systematic error on the extrapolation from the 
simulated 293(5) MeV < f?!j[ <J 417(10) MeV pion mass-range to the physical point. In order to 
address this issue we must simulate with lighter quark masses, which necessitates an increase in 
the physical lattice volume in order to maintain small finite-volume corrections. As increasing 
the number of lattice sites is very costly we must use coarser lattices in order to perform the 
calculation with the currently available resources. Aside from the larger discretization errors, the 
only significant impact of simulating with a coarser lattice is an increase in the size of the residual 
mass m res , which parameterizes the explicit chiral symmetry breaking occuring due to the finite 
length of the fifth dimension. m res gets larger due to the increased number of low-modes of the 
Wilson Dirac operator in the infrared regime, that are likely caused by so-called 'dislocations' - 
localized instanton-like artefacts - in the gauge fields. Configurations containing these low modes 
may be suppressed in the path integral via the introduction to the gauge action of an additional 
weighting factor known as the Dislocation Suppressing Determinant Ratio (DSDR) [|3|-|6j]. 
In this paper we present the '2012 analysis' of the RBC & UKQCD collaboration's /3 = 1.75 
32 3 x 64 x 32 DWF ensembles that make use of the Iwasaki+DSDR gauge action to reach unitary 
pion masses as low as 171(1) MeV and partially-quenched pion masses at a near-physical value 
of 143(1) MeV. The results for the physical quark masses and lattice spacings presented in this 
document were used in our recent calculation of the AI = 3/2 K — > %% amplitudes with physical 
kinematics DVD - 

Note that the pion masses in physical units quoted above and in the abstract, as well as those given 
in the remainder of this paper, were obtained by combining the data at the simulated strange quark 
mass with the final lattice spacings obtained in this analysis, and the error represents the combined 
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systematic and statistical uncertainty. 

Throughout this document we make use of the shorthand 32ID to refer to the 32 3 x 64 x 32 
Iwasaki+DSDR ensemble set, and 321 and 241 for the 32 3 x 64 x 16 and 24 3 x 64 x 16 Iwasaki 
ensemble sets respectively. This notation differs slightly from ref. [[z], where the Iwasaki+DSDR 
ensemble set was labelled 32IDSDR. 

In this paper all dimensionful quantities are expressed in lattice units unless other units are explic- 
itly specified or clarity is served by introducing explicit factors of the lattice spacing a. 
In section |n] we provide further details on the Iwasaki+DSDR gauge action and present the simu- 
lation parameters of our 32ID ensembles. In section [III] we present our results for the pseudoscalar 
masses and decay constants, the Omega baryon mass (used to set the scale), the Sommer scales r$ 
and r\ and also Bk, measured on these ensembles. 

Due to the good chiral symmetry of the domain wall formulation, lattice artefacts involving odd 
powers of the lattice spacing are suppressed and we gain automatic offshell G[d) improvement. 
If we ignore small higher order terms involving products of the square of the lattice spacing and 
the quark mass, the only surviving dependence on the lattice spacing is a single 0(a 2 ) term for 
each measured quantity. Of course this term depends on the lattice action, but as all other parame- 
ters describing the quantity are common between the Iwasaki and Iwasaki+DSDR actions, we can 
easily obtain the a 2 coefficients for the Iwasaki+DSDR action by comparing any single measured 
value on the 32ID ensemble set with the continuum limit obtained from the Iwasaki ensembles. 
In practice we can include the Iwasaki+DSDR ensembles in our simultaneous chiral/continuum 
fitting framework to obtain these a 2 coefficients under a fit to all of the 32ID data, as well as allow- 
ing this data to further constrain the shared parameters. This is discussed further in section|IVl We 
use this procedure in sections IVl through IVIIII to simultaneously fit the aforementioned quantities 
over all three ensemble sets, from which we obtain the lattice spacings and physical quark masses 
as well as improved continuum predictions for the decay constants, Sommer scales and 



II. SIMULATION DETAILS AND ENSEMBLE PROPERTIES 

We generated a set of domain wall fermion ensembles using the Iwasaki+DSDR gauge action, 
which allows for simulations to be performed on coarser lattices while retaining good chiral sym- 
metry and topological tunneling. In this section we provide background on the DSDR term fol- 
lowed by a list of simulation parameters and an analysis of the integrated autocorrelation length 
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and topological charge evolution. 



A. The DSDR Term 



The explicit breaking of chiral symmetry in the domain wall fermion framework can be described 
by an additive mass renormalization parameter referred to as m rcs , whose magnitude is related to 
the eigenvalue density p(A) of the logarithm of the transfer matrix in the fifth-dimension, 

//transfer = 2 tanh 1 f ^ ) , ( 1 ) 

\ 1 + Uw J 

that describes the propagation of quarks through the fifth dimension, via the following relation [8|] : 

m res =/? 4 / dA p(X)e~ LsX . (2) 
Jo 

Here R is a (possibly eigenvalue-dependent) radius factor, Dyy is the Wilson Dirac operator and 
Hw = Y'Dy/ is the hermitian Wilson Dirac operator. 

In the low-eigenvalue region the eigenmodes of //transfer an d those of Hw are necessarily identical. 



It has been demonstrated 1181-11411 that the modes of the latter can be divided into two regions, 
one containing only localized eigenmodes with small eigenvalues and one containing extended 
eigenmodes with large eigenvalues, separated by a mobility edge X c . Picking out the dominant 
contributions above and below the mobility edge from eqn.[2l we expect the following dependence 
of m rcs uponL s : 

m res = R 4 e p ( Xc) +R 4 p(0)^, (3) 

where R e and Rj are the radius parameters for the extended and local modes respectively. The 
exponentially-decreasing contribution from the extended modes above the mobility edge can be 
controlled by increasing L s , with a cost that rises at worst linearly. In our previous Iwasaki simu- 
lations the magnitude of m res was dominated by the term in p(0), the density of near-zero eigen- 
modes. These modes are thought to be associated with localized and short-lived dislocations or 
'tears' in the gauge fields, which can cause changes in the field topology. As the strong coupling 
limit is approached, the gauge fields become more disordered and the density of near-zero modes 
increases sharply. In order to maintain good chiral symmetry properties at stronger coupling we 
must therefore seek to suppress the near-zero modes. On the other hand we must take care not 
to also remove the very-near- zero eigenmodes that are required for topological tunneling to occur 
during the gauge evolution. 
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The DSDR, or 'Auxiliary Determinant' is applied to the gauge action as a multiplicative weight of 
the form yMg] 



W(M;e f ;e b ) 



det [D w ( —M + ie h y 5 ) *D W ( —M + ie h y 5 )} If + ej 



n 



(4) 



det [D w (-M + ie f Y 5 ) f ZV ( -Af + z'£/7 5 )] 
where £y and £^ are tunable parameters with typical sizes < ej <^ £^ < 1. With this weighting, 
the contribution of a single eigenmode to the Molecular Dynamics force becomes a function of £f 
and of the form 



(5) 



which when plotted against the eigenvalue has a peak and tail which are independently tunable 
by varying the two parameters. It is therefore possible to tune the force to suppress near-zero 
eigenmodes while not completely suppressing the essential very-near-zero modes. 
Numerical studies [g] have demonstrated a reduction in chiral symmetry breaking while retaining 
adequate topological tunneling through the use of this term. 



B. Simulation Parameters 



We generated DWF ensembles with the Shamir kernel and the Iwasaki+DSDR gauge action on 
a 32 3 x 64 lattice volume with L s = 32. We used a 'domain wall height' of M5 = 1.8 and a 
gauge coupling of /3 = 1.75, which as determined in section |V] corresponds to an inverse lattice 
spacing of 1.37(1) GeV. We generated two ensembles with bare light-quark masses of m/ = 0.001 
and mi = 0.0042, for which the corresponding unitary pion masses are 171(1) and 246(2) MeV. 
In this document we analyse ~ 1400 and ~ 1200 MD time units on these ensembles respectively 
(discarding 500 and 600 MD time units respectively for thermalization). On each of the ensembles 
we simulated with a single strange-quark mass close to the physical value and use reweighting 
to correct to the true physical value in our fits a posteriori. Further details of the number of 
reweighting steps and stochastic samples are given in the following subsection. 



C. Ensemble Generation 



In this section we provide a summary of the Monte Carlo algorithms that were employed for the 
gauge evolution. Further discussion of our algorithms, along with the full set of parameters, can 
be found in Appendix lAl 
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For the fermionic contribution to the evolution of the m/ = 0.0042 ensemble we employed the 
'RHMC II' algorithm ||15|1 . in which the calculation of the strange-quark determinant is broken 
into three factors and evaluated using the rational approximation with equal molecular dynam- 
ics time steps, and the determinant of the two degenerate light-quarks was preconditioned by the 
strange-quark determinant. With the notation 3>{m) = Dp WF (M5,m)DDWF(^57 m ) for the Her- 
mitian domain wall operator and using M a (m) to represent the rational approximation to the a th 
power of @ for mass m, the algorithm can be written as 



det 



9{m s )^9{mi) 

0(1)2 



det 



9)(m s ) 

~W) 



•det 



•det 



•det 



(6) 



where each determinant is estimated using independent pseudofermion fields. We made use of an 
Omelyan integrator with parameter A =0.22 during the evolution of this ensemble. 
For the lighter m/ = 0.001 ensemble, we were able to achieve a significant speed-up Jji]] in evalu- 
ating the light-quark contribution to the gauge field update using multiple Hasenbusch mass split- 
tings lllq. Il7ll. Here the determinant is split into k steps (with k = 6 in our case), each evaluated 
using a shifted mass: 



det 



k+l 



ndet 



t-1, 



(7) 



where = /io < Mi ■ • -Hk+i = 1 — mi. The intermediate masses /!,•(?' = l..k) can be continuously 
tuned, enabling us to evaluate the individual determinants at a reduced precision - 10 residual as 
opposed to 10~ 8 - considerably reducing the computational cost. The strange-quark determinants 
were again evaluated using the rational approximation. We obtained a further increase in speed by 
utilizing a force gradient integrator [1(3, 18] in place of the Omelyan integrator. 
In table H we give details of the Molecular Dynamics time steps and the update ratios for each 
component of the force, alongside the total MD time, the Metropolis acceptance and the values of 
the average plaquette and chiral condensate on each ensemble. 



D. Ensemble Properties 

In figure CD we plot the Monte Carlo evolution of the plaquette, topological charge and the light- 
quark pseudoscalar density. We measured the topological charge directly using 'cloverleaf ' esti- 
mates of the field strength tensor, with lxl, 1x2, 2x2, 1x3 and 3x3 Wilson loops calculated on 
APE-smeared gauge fields (with 60 smearing steps) and combined using the '51i' (five-loop im- 
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m s 


mi 


m s /rhi 


At x Ar steps 


N G ■ NdSDR : M'erm 


z(MD) 


Acceptance 


(P) 






0.0042 


7.8 


1/8x8 


64:8:(2:1) 


1176 


70% 


0.512198(3) 


0.001579(5) 


0.045 


















0.001 


16.5 


1/9x9 


12:6:1 


1432 


73% 


0.512230(3) 


0.001202(3) 



TABLE I. Simulation parameters for the 32ID ensembles. Here the fifth column contains a gross summary 
of the algorithm, giving the ratio of gauge field updates (No) to the number of DSDR updates (A^dsdr) to 
the number of updates of the fermion force (M' er m)- For the heavier ensemble, the fermion component is 
divided into the rational approximation for the strange-quark determinant and the light quark determinant; 
the former is updated twice as often as the latter. On the lighter ensemble the strange-quark determinant and 
the Hasenbusch-preconditioned light-quark determinant are not nested but instead are evaluated indepen- 
dently and their force contributions combined linearly. The Molecular Dynamics time step for the top-level 
integrator and the number of steps per trajectory (N steps ) is given in the fourth column. The quantity r(MD) 
is the length of the ensemble used for the analyses in this document, measured in Molecular Dynamics time 
units. 



proved) combination Ill9ll which eliminates the £?(a 2 ) and ff(a 4 ) terms at tree-level. We show 
histograms of the topological charge distribution in figure [2l 

Figure |3] contains plots of the integrated autocorrelation time for various quantities on the m\ = 
0.001 and m/ = 0.0042 ensembles as a function of the cut on the upper bound of the integral, A cut : 

Acut 

(8) 



T& tf (A cot ) = |+£C(A) 



A=l 



where 



c ( A) = (im-rxrj+Q-t) 



(9) 



for a quantity Y, where Y is the expectation value over the ensemble, a 2 its variance, and A is the 
Molecular Dynamics time separation between measurements. The average in the second equation 
is performed over the set of pairs of configurations separated by A MD time units. In order to 
correctly estimate the errors on the integrated autocorrelation time, we investigated two strategies: 

1. At each fixed A we formed a bootstrap distribution to estimate the error on the mean (...) t 
in eqn. Prior to performing the bootstrap resampling, we binned the set of measurements 
(Y(t) — Y) (Y(t + A) — Y) over neighboring configurations (indexed here by t). The bin size 
was successively increased until the errors stopped growing, which we found to be at bin 
sizes of 25 and 20 on the m/ = 0.001 and m; = 0.0042 ensembles respectively. The error on 
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500 1000 1500 2000 2500 
MD time units 



500 1000 1500 2000 2500 
MD time units 



FIG. 1. Monte Carlo evolution of the average plaquette (top), topological charge (middle), and light-quark 
pseudoscalar density (bottom) on the mi = 0.001 (left) and m/ = 0.0042 (right) ensembles. 



| □ ml = 0.001 | 



□ ml = 0.0042 



FIG. 2. Topological charge distributions for the m/ = 0.001 (left) and mi = 0.0042 (right) ensembles. 



Tint was obtained from the bootstrap sum over C(A) according to eqn.[8] This method closely 
resembles the standard strategy for binning equivalent quantities over a set of correlated 
measurements under a bootstrap. 



2. We took the full set of measurements Y(t) over the ensemble and formed blocks by aver- 
aging over neighbouring configurations. We then measured the correlations between these 
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blocks, taking the center-point of each block as the associated MD time. This has the ef- 
fect of averaging over short-range correlations, exposing those with longer range, but also 
results in changes to the central value of Tint at fixed A cut as the bin size is increased, as at 
each bin size we are measuring a different quantity. We chose the optimal bin size to be the 
point where further increases resulted in statistically consistent central values. This strategy 
was used in our 2010 analysis for estimating the autocorrelation length of the two Iwasaki 
ensemble sets. 

The aforementioned figure contains plots for both of these strategies. We see that they give con- 
sistent results. The integrated autocorrelation time for the majority of the quantities we looked at 
appears to lie between 5 and 10 MD time units. However, as is typically the case, the topological 
charge (and of course the pseudoscalar condensate) display considerably larger autocorrelation 
lengths, around 25 MD time units on the ligher ensemble and 15 on the heavier ensemble, re- 
flecting their sensitivity to the underlying global gauge field topology. The larger autocorrelation 
length suggests a lower topological tunneling rate for our lighter ensemble. However we empha- 
size that these autocorrelation times are considerably shorter than those of the Iwasaki lattices, 
which were estimated to be ^(80) MD time units yj] from the topological charge measurements. 
For the simulation parameters and properties of the 321 and 241 Iwasaki ensemble sets we refer to 
reader to ref . yj] . 



E. Reweighting the Strange Quark 

We make use of reweighting in the strange sea-quark mass to obtain the mass dependence of our 
data, and hence interpolate to the physical value, without incurring the expense of simulating with 
additional masses. The reweighting factor wi for a particular reweighted mass m™ and configu- 
ration i is determined by measuring the degree to which that configuration, as sampled from the 
un-reweighted path integral, contributes to the path integral with the reweighted mass; in practice 
this involves the calculation of the ratio of Dirac-matrix determinants with the reweighted and 
simulated masses respectively. The expectation value of an observable 6 with the shifted strange- 
quark mass is then obtained by first measuring on the original, unreweighted configurations, then 
applying the reweighting factors: 
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FIG. 3. The integrated autocorrelation time is shown for the average plaquette, topological charge, the chiral 
condensate and pseudoscalar density for the light and heavy quark species (labelled T and 'h' respectively), 
and the pseudoscalar two-point function at t = 20, as a function of the upper bound on the integral A cut , using 
data from the m/ = 0.001 (top) and mi = 0.0042 (bottom) ensembles. For those plots on the left we estimated 
the eiTors by binning the set of correlations between measurements at fixed MD time separation (the first 
strategy discussed in the text), and in those on the right we block over the data and measure the correlation 
between blocks (the second strategy). We chose bin sizes of 25 and 20 on the lighter and heavier ensembles 
respectively. The pseudoscalar two-point function was only measured every 8 MD time units, hence for 
both methods we bin these data with a bin size of 24 MD time units. In the right-hand plots the data have 
been shifted slightly for clarity. 
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The determinants are stochastically evaluated using several Gaussian sampled vectors and the 
weight factor obtained from the average over these samples. This procedure was used in the 2010 
analysis, and more details can be found in ref UJ]. 

We performed measurements over incremental steps from the simulated mass of 0.045 up to 0.052. 
We previously found that the number of stochastic samples required for a reliable estimate of 
the weighting factor is dependent upon the size of the mass increments, with smaller increments 
requiring less samples. As a result, we use two stochastic samples and small increments of Am/, = 
0.00025 - the same parameters as were used for the 241 ensembles. 

The reweighting procedure naturally reduces the effective number of configurations 7V e ff in each 
ensemble set. In ref yj] we showed that a reliable estimate of this quantity can be determined via 
the following expression: 



A value of unity indicates that the measurement is entirely dominated by a single configuration, 
whereas N e ff is equal to the original number of configurations Af con f when there are no fluctuation 
in the weighting factors. In section |V] we measure the physical strange quark mass to be m^ hys = 
0.0467(6), which is close to the simulated value. At the nearest reweighted mass-step to the 
physical mass, that with m h = 0.0465, we find N e g = 133 (A^onf = 180) and N e ff =119 (N con f = 
148) on the /«/ = 0.001 and /«/ = 0.0042 ensembles respectively, suggesting that reweighting to 
the physical strange-quark mass will result in only a 10-15% increase in the statistical errors on 
these ensembles. This is of a similar magnitude to the increase suggested by the values of 7Y e ff 
on the 321 ensembles, which are given in ref Jl]. On the 241 ensembles we require a slightly 
larger extrapolation to reach the physical value, hence the reweighting introduces larger increases 
of 25-35%. 



HI. RESULTS FROM THE 32 3 DWF+ID ENSEMBLES 



In this section we present the results of fitting to a number of observables on the 32ID ensembles. 
We performed measurements on 180 configurations on the m/ = 0.001 ensemble and 148 on the 
mi = 0.0042 ensemble, with each configuration separated by 8 MD time units. The analysis in the 
previous section suggests an autocorrelation length of ~ 25 on the m/ = 0.001 ensemble and ~ 7 
on the mi = 0.0042 ensemble, which can be overcome by binning the data before performing the 
fits. We shifted the gauge fields in the time-direction by 16 lattice spacings relative to the previous 
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configuration prior to measuring the quark propagators. This has the effect of reducing the corre- 
lation between successive measurements, suggesting that binning the data may not be necessary. 
However, this does not apply to the measurements of the Sommer scales r$ and r\, which are 
formed using Wilson loops with origins on all lattice sites. In order to remain consistent, we de- 
cided to bin the data for all of our quantities over 4 successive measurements (32 MD time units) 
on both ensembles; although this is larger than the measured autocorrelation length, it matches 
the periodicity of the quark propagator measurements, and is therefore a more natural choice. We 
found no statistically significant dependence on the bin size in any of our measured error values, 
hence the choice of bin size has little effect on the final results of this analysis. 
The pseudoscalar meson two-point correlation functions were calculated in the same manner as 
those on the 321 ensembles, namely using Coulomb gauge-fixed wall source propagators origi- 
nating at the lattice time boundary t = with both periodic (p) and anti-periodic (a) boundary 
conditions in the temporal direction. Taking the p + a combination of propagators to form each 
leg of the correlation function projects out the component travelling forwards in time. Likewise, 
the p — a combination projects out the degenerate backwards-propagating state. The correlation 
functions formed using these combinations of propagators have a temporal periodicity of double 
the usual length, which results in a significant reduction in round-the-world propagation. The 
Omega baryon correlation functions were calculated separately using box-sources with a spatial 
volume of 15 3 lattice sites and with one corner at the spatial origin. These were placed on time- 
slices t = and 32, and anti-periodic boundary conditions were used for the propagators. As 
mentioned above, the gauge fields were shifted in time by 16 units with respect to the previous 
configuration prior to performing all of these measurements. 

For each quantity we tabulate the results of fitting to the time-dependence of the corresponding 
correlation functions measured at the simulated strange-quark mass, and we present example ef- 
fective mass plots demonstrating the quality of our data. We also provide tables of data corrected 
to the physical strange-quark mass of m s = 0.0467(6) determined in section[Vj using the the NLO 
ChPT with finite-volume corrections parameterization for the mass dependence. 
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m x 


mi 




0.001 


0.0042 


0.0001 


0.0018447(60) 


0.0018888(48) 


0.001 


0.0018510(43) 


0.0018889(47) 


0.0042 


0.0018269(58) 


0.0018735(48) 


0.008 


0.0018025(57) 


0.0018500(48) 


0.035 


0.0016939(44) 


0.0017356(39) 


0.045 


0.0016739(39) 


0.0017141(37) 


0.055 


0.0016619(36) 


0.0017014(35) 



TABLE II. m res on the 32ID ensemble set at the simulated strange-quark mass. 
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FIG. 4. The chiral extrapolation of m! Tes over the unitary data points at the simulated strange-quark mass, 
with m res in the chiral limit denoted by the brown square point (left); and the strange-quark mass dependence 
of m res in the chiral limit (right). 

1. The Residual Mass 



The residual mass at a non-zero (partially-quenched) quark mass m x may be determined via the 
following ratio: 



(0|J?|w) 



(12) 



where /| is the pseudoscalar density at the midpoint of the fifth dimension, and 7| is the physical 
pseudoscalar density constructed from the surface fields. The prime superscript is used to differ- 
entiate this quantity from the residual mass in the two-flavor chiral limit, m res = m[Jm x = m/ = 0). 
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We averaged the data at t and T — t (we refer to this as folding the data), where T is the lattice 
temporal extent, and fit over the time range 6-32 on both ensembles, obtaining the values given in 
table HIl Note that on the lighter ensemble, the non-unitary values were determined on a reduced 
data set of 92 configurations; these data were not used in the later analysis but are presented here 
for completeness. We obtained m lcs by extrapolating over the unitary light-quark mass to the chiral 
limit at each reweighted strange-quark mass. A plot of the chiral extrapolation at the simulated 
strange-quark mass is shown in figure HI Owing to the minor strange-quark mass dependence of 
this quantity evident in the right panel of figure HI and the small separation between the simu- 
lated and physical strange quark masses, the value of m rcs at the physical strange quark mass is not 
measurably different from that at the simulated value of 0.001842(7). 



2. Pseudoscalar Masses 



We calculated a series of pseudoscalar meson two-point functions of the form: 

«S^(0 = <0|O?(0O?(0)|0>. (13) 

Here the subscripts index the interpolating operators and the superscripts denote the operator 
smearing (wall W or local L) at the sink and source respectively. In the following we refer to 
these by the shorthand OiO^ 2 , for example using AA LW to denote the axial-axial correlator with 
wall source and point sink. The pseudoscalar masses were determined via a combined fit to the 
following five correlation functions: PP LW , AP LW and AA m , PP WW and AP WL . The correlation 
functions exhibit the following time dependence: 



°i°* w 2 mxy V 



(14) 



where the sign in the square brackets is + for the PP and AA correlators and — for the AP correla- 
tors. We denote the amplitudes as 



Taking full advantage of the doubled time-extent of the lattice, we performed our fits over the 
time range 8-63 on both ensembles, obtaining the masses listed in table JVJ The values at the 
(unitary) physical strange-quark mass are given in tables [V] and [VI] for the light-light (pion-like) 
and strange-light (kaon-like) quark mass combinations respectively. In figures [5] and [6] we show 
example effective mass plots for the data at the simulated strange-quark mass on the m/ = 0.001 
ensemble. 
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3. Pseudoscalar Decay Constants 

The pseudoscalar decay constants fxy were calculated from the two-point function amplitudes via 
the following equation: 

/ 2 ,A^ Lwl 

s ^H^d^- (16) 

Here Za relates the local four-dimensional axial current - formed with the domain wall surface 
fields - to the Symanzik-improved axial current A s £, and thus renormalizes the local current into 
the continuum normalization. The indices a and /i correspond to the flavor and Euclidean direction 
respectively. 

For domain wall fermions, a partially-conserved five-dimensional axial current can also be de- 
fined, which is related to the Symanzik improved current by a different renormalization coefficient 
Zjrf. Prior to the 2010 analysis, it was typically assumed that the difference between Z^ and unity 
was negligable, hence Za was assumed equal to Za/Z^. This can be obtained using the improved 
ratio Il20ll of the partially-conserved 5D axial current matrix element (j^(f)P(O)) to the local axial 
current matrix element (A$(t)P(Qi)) . As discussed in ref. 112 111 , the assumption that Z^ = 1 is only 
true up to terms &{am )\, leading to an additional &{!%) systematic error in our earlier results. 



However, in refs. 12211 and (12111 it was shown that Za is approximately equal to Zy - the ratio 

to the local vector current V£ 



of the Symanzik-improved vector current V$ a to the local vector current Vf[ - with their differ- 



ence of order mjr es . Since the ratio Zy of the conserved 5D domain wall vector current to its 
Symanzik-improved current is unity up to terms &{a 2 ), this led to the observation that Za 
can be determined much more accurately via the ratio of the local and 5D vector currents, Zy /Zy, 
calculated using the following expression: 

z v _ i? = iigTTOOv; fl (5,o) 
z v LU^vt(x,t)vt(o,o) 

in the limit t 3> a. We calculated Zy jZ-y on 192 and 93 configurations of the m\ = 0.001 and 
0.0042 ensembles respectively, and fit to folded data over the time intervals 8-12 and 7-17. Fig- 
ure U\ shows Zy /Zy(m x = ni[ = 0.001) as a function of time, illustrating the quality of our data. In 
the same figure we also show the chiral extrapolation of the results to m/ = — m res - In table Hill we 
give the fit results on both ensembles and the chirally extrapolated values. For completeness we 
also calculate the ratio Za/Zj^ using the aforementioned ratio [2(1, fitting over the time interval 
5-30 to folded data. The values of this quantity on each ensemble and in the chiral limit are also 
given in table Unl and we show an example correlation function in figure [7] alongside a plot of 
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m q Z A /Z^ Z v /Zy 

0.0042 0.68901(9) 0.6637(46) 
0.001 0.68828(15) 0.6685(36) 
-m res 0.68778(34) 0.6728(80) 

TABLE III. Results for Z^/Z^ and Zy/Zy at the simulated strange-quark mass. 



the chiral extrapolation to mj = — m res . The value of Z^/Z^ at the physical strange-quark mass is 
indistinguishable from the value at the simulated mass. Currently we have not measured Zy/Zy 
on reweighted configurations, however the lack of measurable strange-quark mass dependence of 
Za/Z^ suggests this will not have any effect on our conclusions. 

We calculated the normalized decay constants using the above ratios. The values at the simulated 
strange-quark mass are listed in the second column of table [IV] and the pion-like and kaon-like 
decay constants at the physical strange-quark mass are given in the second columns of tables [V] 
and [VI] respectively. For these quantities, the statistical uncertainty on Zy/Zf is considerably 
larger than that of the bare decay constant. For example, the bare unitary value on the m/ = 0.001 
ensemble has a 0.3% error compared to 1.2% on the normalized quantity. The error on Z^jZ^ 
is much smaller, and if used to normalise the bare decay constants has virtually no effect on the 
relative error. However we chose to continue using Zy/Zy to normalize the decay constants in 
order to eliminate the systematic error associated with using the axial currents. 



4. Omega Baryon Mass 

We determined the Omega baryon masses using box-source propagators with antiperiodic bound- 
ary conditions. In order to improve our statistics we averaged the degenerate upper and lower 
spin-components of the correlation functions prior to fitting. Our fits were performed over the 
interval t = 3-10 on both ensembles, giving the values listed in table I VIII In figure[8]we show the 
effective mass of the Omega baryon on the m/ = 0.001 ensemble with w/, = m x = 0.045, demon- 
strating the quality of our data. 
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5. Neutral Kaon Mixing Parameter 

The neutral-kaon mixing parameter was obtained by fitting the time dependence of the follow- 
ing correlation function to a constant: 

g lat m _ (K°(h)\0 WW+AA (t)\K O (t 2 )) 

K[> %(K°{t l )\A (t))<Ao{t)\K<>(t2)y 

where i^Vv+aa is the AS = 2 four-quark operator responsible for the mixing. This operator is 
inserted at all times t between t\ and t 2 . We form the forwards-propagating K° state using the p + a 
combination of propagators, and the backwards-propagating K° state using the p — a combination; 
in effect this sets t\ = and t 2 = 64 and reduces the round-the-world effects associated with the 
kaons propagating through the temporal boundaries. We performed our fits over the time interval 
8-56, giving the values listed in table I Villi We show an example matrix element in figure[8]and list 
the values of Bxy at the physical strange-quark mass in table|IXl Note that Bk is a renormalization- 
scheme dependent quantity and must therefore be renormalized into a common scheme prior to 
being included in our simultaneous fits; this is discussed in more detail in section [VTT1 

6. The Sommer Scales 

Finally, we obtain the Sommer scales ro and r\ using Wilson loops formed from products of 
time-directed gauge links, for which closure is not required due to Coloumb gauge-fixing. The 
time dependence of the Wilson loop W(r,t) was fit from t = 3 to 8 for each value of the spatial 
separation r, and the resulting potential V(r) then fit over the range r = 2.00 — 9 to the Cornell 



potential D23Q 



cc 

V(r)=V -- + ar, (19) 
where Vq, a and o are constants. The Sommer scales are determined directly from the potential: 



where Aq = 1.65 and A\ = 1.00 for ro and r\ respectively. In figure [9] we show an example of 
the effective potential V e s(t) at r = 2.45 on the m/ = 0.001 ensemble and the resulting fit to the 
potential V(r) using the Cornell form. In table |X] we give the values of r\ and ro, as well as their 
ratios, at the simulated and physical strange-quark masses. 
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m x 


My 


(0.001) 


m^, (0.0042) 


U (o.ooi) 


/^(Q.0042) 


0.055 


0.055 


0.5463(2) 


0.5476(2) 


0.1354< 


16) 


0.1363( 


16) 


0.045 


0.055 


0.5207(2) 


0.5220(2) 


0.1324( 


16) 


0.1334( 


16) 


0.035 


0.055 


0.4941(2) 


0.4954(2) 


0.129K 


16) 


0.1302< 


16) 


0.008 


0.055 


0.4159(3) 


0.4174(3) 


0.1183( 


14) 


0.1200( 


15) 


0.0042 


0.055 


0.4041(4) 


0.4055(4) 


0.1164< 


14) 


0.1184( 


15) 


0.001 


0.055 


0.3942(6) 


0.3955(6) 


0.115K 


14) 


0.1173( 


15) 


0.0001 


0.055 


0.3915(7) 


0.3928(7) 


0.1149( 


14) 


0.1173( 


15) 


0.045 


0.045 


0.4940(2) 


0.4953(2) 


0.1294( 


16) 


0.13051 


16) 


0.035 


0.045 


0.4662(2) 


0.4675(2) 


0.1262( 


15) 


0.1274( 


15) 


0.008 


0.045 


0.3831(3) 


0.3846(3) 


0.11561 


14) 


0.1174( 


14) 


0.0042 


0.045 


0.3703(3) 


0.3718(4) 


0.1137( 


14) 


0.1158( 


14) 


0.001 


0.045 


0.3594(5) 


0.3610(5) 


0.1123( 


14) 


0.1148( 


14) 


0.0001 


0.045 


0.3564(6) 


0.3581(6) 


0.112K 


14) 


0.1147( 


15) 


0.035 


0.035 


0.4368(2) 


0.4381(2) 


0.123K 


15) 


0.1243( 


15) 


0.008 


0.035 


0.3476(3) 


0.3491(3) 


0.11261 


14) 


0.1144( 


14) 


0.0042 


0.035 


0.3334(3) 


0.3350(3) 


0.11071 


13) 


0.1129( 


14) 


0.001 


0.035 


0.3212(4) 


0.3230(4) 


0.1092( 


13) 


0.1118( 


14) 


0.0001 


0.035 


0.3178(5) 


0.3197(5) 


0.1090( 


13) 


0.1118( 


14) 


0.008 


0.008 


0.2273(2) 


0.2287(3) 


0.1024( 


12) 


0.1044( 


13) 


0.0042 


0.008 


0.2048(2) 


0.2063(3) 


0.1005( 


12) 


0.1028( 


13) 


0.001 


0.008 


0.1839(2) 


0.1854(3) 


0.09881 


12) 


0.1015( 


12) 


0.0001 


0.008 


0.1775(2) 


0.1791(3) 


0.0984( 


12) 


0.1013( 


13) 


0.0042 0.0042 


0.1795(2) 


0.1810(2) 


0.0986( 


12) 


0.101K 


12) 


0.001 


0.0042 


0.1549(2) 


0.1564(2) 


0.0969( 


12) 


0.0997( 


12) 


0.0001 0.0042 


0.1472(2) 


0.1487(3) 


0.0964( 


12) 


0.0994) 


12) 


0.001 


0.001 


0.1250(2) 


0.1265(2) 


0.09501 


12) 


0.098 1( 


12) 


0.0001 


0.001 


0.1151(2) 


0.1167(3) 


0.0944( 


12) 


0.09771 


12) 


0.0001 0.0001 


0.1042(2) 


0.1058(3) 


0.0938( 


12) 


0.0973( 


12) 



TABLE IV. Pseudoscalar masses m^ (mi) and decay constants f xy (mi) on the 32ID ensembles at the simu- 
lated strange-quark mass (m/, = 0.045). 
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m x 


m y 


rn^y (0.001) 


ntxy (0.0042) 


^,(0.001) 4(0.0042) 


0.008 


0.008 


0.2272(2) 


0.2286(3) 


0.1027(12) 0.1048(13) 


0.0042 


0.008 


0.2048(2) 


0.2063(3) 


0.1008(12) 0.1031(13) 


0.001 


0.008 


0.1838(2) 


0.1854(3) 


0.0991(12) 0.1018(13) 


0.0001 


0.008 


0.1775(2) 


0.1791(3) 


0.0987(12) 0.1016(13) 


0.0042 0.0042 


0.1794(2) 


0.1809(3) 


0.0989(12) 0.1014(12) 


0.001 


0.0042 


0.1548(2) 


0.1563(3) 


0.0972(12) 0.1000(12) 


0.0001 0.0042 


0.1471(2) 


0.1486(3) 


0.0967(12) 0.0997(12) 


0.001 


0.001 


0.1249(2) 


0.1265(3) 


0.0953(12) 0.0984(12) 


0.0001 


0.001 


0.1151(2) 


0.1166(3) 


0.0947(12) 0.0981(12) 


0.0001 0.0001 


0.1042(2) 


0.1058(3) 


0.0941(12) 0.0976(12) 



TABLE V. Pion masses m^ (mi) and decay constants f xy (mi) on the 32ID ensembles at the physical strange- 
quark mass (nth = 0.0467(6)). 



m x 


mj* (0.001) 


m xh (0.0042) 


fxh (0.001) f xh (0.0042) 


0.008 


0.3890(20) 


0.3903(21) 


0.1161(14) 0.1178(15) 


0.0042 


0.3762(21) 


0.3777(22) 


0.1143(14) 0.1163(15) 


0.001 


0.3653(21) 


0.3669(23) 


0.1128(14) 0.1153(15) 


0.0001 


0.3624(21) 


0.3640(23) 


0.1126(14) 0.1153(15) 



TABLE VI. Kaon masses m x h (mi) and decay constants f x h (mi) on the 32ID ensembles at the physical 
strange-quark mass (m/, = 0.0467(6)). 





m h 


m a (0.001) 


m a (0.0042) 


0.055 


0.045 


1.2641(34) 


1.2735(36) 


0.045 


0.045 


1.2130(37) 


1.2220(41) 


0.035 


0.045 


1.1608(42) 


1.1695(48) 


0.0467(6) 0.0467(6) 


1.2248(77) 


1.2326(55) 



TABLE VII. Omega baryon masses on the 32ID ensembles at the simulated strange quark mass mj, = 0.045 
(first three rows) and at the physical strange quark mass (fourth row). 
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m x 


nty 


B r) .(0.001) 5^(0.0042) 


0.008 


0.055 


0.645(2) 


0.645(2) 


0.0042 


0.055 


0.643(4) 


0.645(4) 


0.001 


0.055 


0.650(16) 


0.665(16) 


0.0001 


0.055 


0.665(28) 


0.689(28) 


0.008 


0.045 


0.629(2) 


0.628(1) 


0.0042 


0.045 


0.626(3) 


0.625(3) 


0.001 


0.045 


0.630(10) 


0.632(10) 


0.0001 


0.045 


0.639(17) 


0.644(18) 


0.008 


0.035 


0.610(1) 


0.609(1) 


0.0042 


0.035 


0.605(2) 


0.604(2) 


0.001 


0.035 


0.606(6) 


0.602(6) 


0.0001 


0.035 


0.610(10) 


0.605(10) 



TABLE VIII. The partially-quenched neutral kaon mixing parameter (mi) on the 32ID ensembles at the 
simulated strange-quark mass (m/, = 0.045). 



m x 


B xh (0.001) 


5^(0.0042) 


0.008 


0.632(2) 


0.631(2) 


0.0042 


0.630(4) 


0.628(3) 


0.001 


0.638(11) 


0.635(11) 


0.0001 


0.651(17) 


0.649(19) 



TABLE IX. The partially-quenched neutral kaon mixing parameter m x k (mi) on the 32ID ensembles at the 
physical strange-quark mass (nth = 0.0467(6)). 
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mi 


ro r x ri/r 


0.045 

0.0042 

0.0467(6) 


3.2732(63) 2.1208(97) 0.6479(36) 
3.2616(75) 2.1270(105) 0.6521(37) 


0.045 

0.001 

0.0467(6) 


3.2977(62) 2.1346(98) 0.6473(34) 
3.2959(73) 2.1401(100) 0.6493(37) 



TABLE X. The Sommer scales ro and r\ and their ratio on the 32ID ensembles at the simulated strange 
quark mass nif, = 0.045 (first and third rows) and at the physical strange quark mass (second and fourth 
row). 
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FIG. 5. Effective unitary pion masses on the mi = 0.001 ensemble from the PP LW correlator (top left), 
PP WW correlator (top right), AP LW correlator (center left), AP WW (center right) and AA LW correlator 
(bottom). Note the different vertical scale for the WW correlators. The horizontal bands represent the result 
for the mass from a simultaneous fit. 
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FIG. 6. Effective unitary kaon masses on the m/ = 0.001 ensemble from the PP LW correlator (top left), 
PP WW correlator (top right), AP LW correlator (center left), AP WW (center right) and AA LW correlator 
(bottom). Note the different vertical scale for the WW correlators. The horizontal bands represent the result 
for the mass from a simultaneous fit. 
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FIG. 7. Zy/Zy (top left) and Z^/Z^ (top right) as a function of time, calculated with unitary quarks on 
the mi = 0.001 ensemble. The bottom figure shows the chiral extrapolation of Zy /Zy and Za/Z^. In these 
plots the ratios have been abbreviated to Zy and Z4. 




12 3 4 



FIG. 8. The left panel displays the fit to the Q baryon mass with valence strange mass m x = 0.045 on the 
mi = 0.001, m/j = 0.045 ensemble on the 32ID lattice, showing the quality of the fit with our box source. 
The right panel shows the B^, matrix element with m x = m y = 0.001 as a function of time on the same 
ensemble. 
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FIG. 9. The left panel shows the effective potential of the Wilson loops with a spatial extent of r = 2.45 on 
the mi = 0.001 ensemble at the simulated strange-quark mass, overlaid by the fit to the range t = 3-8. The 
right panel shows the static inter-quark potential V (r) on this ensemble, again at the simulated strange-quark 
mass, as a function of the spatial extent of the Wilson loops, overlaid by the fit to the Cornell form over the 
range r = 2.00-9.00. 
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IV. SIMULTANEOUS CHIRAL/CONTINUUM FITTING PROCEDURE 

In order to extrapolate to the continuum limit and physical quark masses we perform a simultane- 
ous global fit over our three ensemble sets. In this section we detail the fitting procedure and the 
subsequent chiral/continuum extrapolation. In addition we discuss the differences between this 
analysis and the 2010 analysis Jl],Q] of the 24 3 and 32 3 DWF+I ensembles. 



A. Global Fits and Scaling 

For a given choice of lattice action and a given bare coupling /3 , 2 + 1 flavor lattice QCD has two 
free parameters: the relevant couplings representing the quark masses. For 2+1 flavor QCD these 
are the average up/down quark mass am u / d and the strange quark mass m s , expressed in lattice 
units. We can picture taking the continuum limit of the discretized theory as gradually taking 
/3 — > oo while following a curve of am u / d (fi) and am s (f5) that fixes the continuum physics to that 
of the real world; this curve is known as a scaling trajectory. Experimental inputs are used to 
determine the lattice spacing and physical quark masses for each bare coupling, and this imposes 
a constraint on each point on this scaling trajectory. (Our standard choice is to require that m^, 
m %/ m Q. and mx/niQ take their physical values.) This in turn allows us to constrain the continuum 
limit we determine to be the physical point. 

We can relate two points (am/,am/,,/3) and (a'm'^a'm'^fi') that lie on a particular scaling trajectory 
via two scaling parameters Z\ and Z/ 2 , defined as yj]: 



z ' W)= Wll <21) 



where / E {l,h}. Here 



W) = ^ (22) 

is the ratio of lattice spacings and rhf = mf + m KS , where m res is the residual mass of domain 
wall QCD. In practice, we define our scaling parameters using the /3 = 2.25 (321) ensemble as 
a reference; we refer to this as the primary ensemble set, on which Z/, Z/, and R a are unity by 
definition. We may interpret our matching of quark masses to the bare masses on our primary 
ensemble set as a convenient, if inelegant, intermediate renormalization scheme, for which the 
regularization involves an explicit choice of lattice action and bare coupling, and whose values are 
determined by the hadronic inputs. The renormalization scale in this scheme is the scale at which 



28 



the bare mass is defined: the inverse lattice spacing of the primary ensemble. The renormalized 
masses are then 

m r { =Z l R a [a'm' l }/a, (23) 
m r l = Z h R a [a'm' h ] /a , (24) 

where unprimed quantities are defined on the primary ensemble set and primed quantities on some 
other ensemble set, and a' is related to a via a' = R a [ a. 

Considering only the unitary observables for simplicity, any observable Q is a function of the bare 
quark masses and the bare coupling. We take this as Q(a > in' I1 a'm' h , /3 ; ), at coupling and quark 
masses differing from the primary ensemble set. This can equally be expressed as a function of 
the renormalized quark masses and the lattice spacing as 

Q(a'm h dm h ,P') = /(mf, m£, a 2 ) . (25) 

The function / depends on the lattice action and on the choice of physical quantities used to 
determine the scaling trajectory. Since among these input parameters is a quantity with a physical 
scale (in our case the Qr mass), we choose to view the function as depending on this scale so its 
arguments can be expressed in physical units. The function itself will have a continuum limit as /3 
and ft' become large. 

Consider a double expansion in quark masses and in lattice spacing around our primary ensemble 

/(mf , m£', a' 2 ) = f(m h m hl a 2 ) 

+ ^«-%) (26) 
+ ^(a' 2 -a 2 ). 
+ 0(m 2 ,a 2 fhf,a A ) 

where the partial derivatives are evaluated at R a = Z/ = Z/, = 1 . 

If / is a quantity used to determine the scaling trajectory then we necessarily constrain that j^t = 
at the match point. In this paper we introduce a new DSDR term to the effective gauge action. To 
this order only the term j^j depends on the the lattice action. We can therefore determine the 
parameters of / for a given parameterization, accurate to this order, via a fit to a set of points 
over multiple ensembles, and including the two different gauge actions. Even though there is only 
a single lattice spacing with the DSDR gauge action, it will usefully contribute to a universality 
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constrained global fit by significantly constraining the mass dependent terms in a global parame- 
terization of /(m/, m/ 7 , a 2 ). 

For the purposes of matching ensemble sets with different lattice spacing we ignore terms of higher 
order in dmj, 5m/, and in a 2 . Since we allow Z/, for masses in the region of the strange quark to 
differ from Z/ for masses in the region of the up/down quarks, in this matching context we may 
consider small variations in the quark masses only. 

We can also see immediately that if we instead use a nearby reference point (am/, am/,, j8) — > 
(a[mi + A;],a[m/ 2 + A/,], /3), the ratios Zf and R a change only by terms 6 («(j8) 2 — a'(p') 2 ) with 
coefficients that are functions of the mass differences Ay that vanish as Af — > 0. This means that 
there is an allowed range overwhich Z/ and Z/ 2 may be simply taken as a constant. Higher order 
terms in quark masses are, of course, subsequently introduced in our global chiral-continuum fits, 
and we introduce Z^ and Z/ as free fit parameters multiplying quark masses in the allowed range. 
In practice we even find Z/ ~ Z/,; were the matching and primary ensemble sets taken sufficiently 
close to the continuum limit, such that lattice artifacts were small and 1/a, then we would 

necessarily find Z/ = Z/ 2 . The matching scheme can therefore be considered mass -independent 
as the mass dependence of the renormalization factors drops out when the renormalization scale 
becomes large. 

In the following subsection we discuss our strategy for determining the scaling parameters Z;, Z/ 7 
and R a . 



B. Determination of the Scaling Parameters 

In our analysis yj] of the 24 3 and 32 3 DWF+I ensembles, we determined R a , Zi and Z^ by matching 
our lattice data at an unphysical light and heavy quark mass within the range of available data 
on the two simulations. The matching was performed by first choosing a suitable match point 
on one of the ensemble sets (labelled M) which can, but does not necessarily have to be, the 
primary ensemble. On every other ensemble set e (in the 2010 analysis only the 241 ensemble set 
remained), two dimensionless ratios, Ri = mu/nihhh and Rh = mih/m^h, were linearly interpolated 
in the unitary light and heavy quark masses until their values matched those measured at the match 
point on ensemble M. Here mu, m\] x and m^hh are respectively the pion, kaon and Omega baryon 
masses measured at unphysical light (Z) and strange (h) quark masses. The match point was chosen 
to minimise the distance of interpolation required on the ensemble sets e. This procedure provides 
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a pair of equivalent masses (in lattice units), (am/) e and (am/,) e , for each ensemble set. Using 
these masses we determined Z/ and Z/, using eqn. [2T1 calculating R a from the ratio of the Omega 
baryon masses at the match point: 

z e = 1 ifl^f) 1 R* = — = m hhh( M J^h) (27) 

where / G {/,/*}. 

The above procedure defines the scaling parameters such that ra//, m Ux and scale perfectly up 
to terms ff{ma 2 ) within the allowed region around the match point. Note that this choice is not 
unique; we could for instance use the pion and kaon decay constants, /// and fa, and the Sommer 
scale ro, and match ro/// and ro///j. R a can then be determined from the ratio of ro measured at the 
match point on each ensemble set. In this case ro, /// and fu x would have no ff(a 2 ) dependence 
instead. In ref. [1] we demonstrated that this produces results that are completely consistent. 
The benefit of this fixed trajectory method is that it enables the separation of the matching from the 
complexities of the subsequent global fits. However, in our combined analysis of the DWF+I and 
DWF+ID ensemble sets, we find that, apart from the lightest partially-quenched point, the range 
of light quark masses on the 241 ensemble set does not overlap with that on the 32ID ensemble 
set (cf. figure [H]). As a result, matching the 241 and 32ID ensemble sets to the 321 primary 
ensemble set at a single point would require a long extrapolation beyond the unitary mass range. 
In addition, the use of independent linear interpolations on each ensemble set is more vulnerable 
to statistical fluctuations than if we were to fit over all data simultaneously. As a result we choose 
the alternate generic scaling method [1], in which R a , Z/ and Z/, are left as free parameters which 
are determined, alongside the low-energy constants, in a global fit of m n , niK and ma over all 
ensemble sets. Here the three conditions that define the scaling trajectory are imposed by omitting 
scaling terms up to &{ma 2 ) from the fit forms describing these quantities, and the values of the 
ratios are selected as those that minimize the global % 2 . In ref. yj] we demonstrated that this 
approach gives consistent results with the fixed trajectory approach. 

Prior to discussing our fit ansatze, it is illustrative to compare the ratios of various dimensionless 
quantities between the 321 and 32ID ensemble sets at a particular match point, using the scal- 
ing parameters determined later in section 13 This allows us to visualize the magnitude of the 
scaling corrections for each quantity. Choosing [am;] 321 = 0.004 and the physical strange-quark 
mass [arrih] 321 = 0.0263(10) as a match point, we used the scaling parameters listed in table IXIIIl 
(combining the statistical error with the systematic errors determined using the procedure given in 
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FIG. 10. Ratios of various dimensionless combinations of observables between the 321 and 32ID ensemble 
sets. The combination of physical quantities is given above or below the corresponding point. A ratio of 
unity indicates perfect scaling between the two ensemble sets. 



section |VBl) to determine the corresponding point on the 32ID ensemble as [am/] 3210 = 0.0066(3) 
and [am/j] = 0.0467(6). We then performed linear fits to a range of quantities over each ensem- 
ble set independently and interpolated each to the corresponding match point quark masses. In 
figure [10] we plot the ratio of a number of dimensionless combinations of these quantities between 
the two ensemble sets. It is immediately apparent that the scaling parameters do indeed fix m K , 
niK and to scale between the two ensemble sets, and the errors on the ratios of these quantities 
are indicative of the size of higher order corrections - in a fixed trajectory matching at this point 
those errors would be zero by definition. 

Considering combinations of m K , m# and with other quantities that retain a scale dependence, 
and for the purpose of making a crude estimate ignoring the small discretization error on the 321 
measurements, we can use this plot to read off the rough size of the discretization error for the 
measurements on the 32ID ensemble set: we estimate 0(3-5%) discretization terms for fu and 
fih, 0(1-2%) for ro, and then a slightly larger 0(5-1%) contribution for r\. 
As an aid to the reader, we also use the aforementioned scaling parameters to place all of the 
simulated quark masses on a common scale, and draw a line to indicate the physical point as 
determined in section |V] These plots are shown in figure ITT1 
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FIG. 1 1. Simulated quark masses on each of our three ensemble sets brought into a common normalization 
with the bare quark masses on our 321 ensemble set using the scaling factors determined in section El The 
top panel shows the light quark mass regime and the bottom panel the heavy quark mass regime. Circular 
points are used to mark the unitary masses and square points the partially-quenched masses. The physical 
up/down and strange quark masses are marked with dashed lines. 



C. Chiral/Continuum Fitting Strategy 

The chiral/continuum fit forms are obtained via a joint expansion in a 2 and ifif. As in ref. 
we consider both an NLO expansion around the SU(2) chiral limit using partially-quenched chiral 
perturbation theory (PQChPT) and also a leading-order analytic expansion about an unphysical 
light-quark mass. Including finite- volume effects in the ChPT, this provides three fit ansatze, which 
we label 'analytic', 'ChPT' and 'ChPTFV, where the latter two refer to the chiral perturbation 
theory forms without and with finite-volume corrections respectively. For each ansatz we expand 
the heavy-quark mass dependence to linear order in the vicinity of the physical strange-quark mass. 
We use a power-counting scheme whereby terms of order irifa 2 and higher are neglected. This 
truncation leaves only a single a 2 term arising from the expansion of the leading order parameter. 
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For example, the analytic form for the pion decay constant /// in physical units is as follows: 

/„ = C{? ( 1 + Cta 2 ) + C{* (m* - mf ) + (mf - mf ) + C 3 /jt (mf - mf ) , (28) 

where the superscript R indicates a renormalized physical quark mass (in a general scheme), and 
mf Q and m^ are the expansion points for the light and heavy quark masses respectively. In our 
power counting scheme, a term in the lattice spacing arises only in the expansion of the leading 
term Cq . It is important to note that the a 2 coefficients parameterizing the lattice artifacts will 
differ between the Iwasaki and Iwasaki+DSDR gauge actions, therefore for the remainder of this 
work we label these coefficients with a superscript denoting the lattice action. 
As discussed in ref. yj], the scaling parameters Zf and Zf that relate the quark masses between 
the ensemble e and the primary ensemble set can be thought of as 'renormalization coefficients', 
removing the ultraviolet divergence and converting the masses into a mass-independent 'matching 
scheme' defined with lattice regularization at /3 = 2.25. It is therefore unnecessary to renormalize 
the input quark masses into a continuum renormalization scheme such as ms prior to performing 
the fits; we need only convert the input masses into the matching scheme. The predictions for the 
physical up/down and strange quark masses can be converted into a more conventional scheme a 
posteriori; this is performed in section IVTl 

In the 2010 analysis, we performed our fits to quantities in physical units. However this required 
us to continually update the lattice spacings and physical quark masses based on the results of 
the fit, iterating until convergence. For this analysis we instead fit to quantities in lattice units, 
which removes the need to repeat the global fit multiple times. However, for clarity, we continue 
to quote our fit forms in dimensionful units; the correctly normalized versions in lattice units for 
an ensemble e can easily be obtained by inserting powers of a e where appropriate to make the 
measurement and input quark masses dimensionless; applying factors of Zf and Zf as before to 
bring the quark masses into the normalization of the primary ensemble; substituting a e with a 1 /R%; 
and finally setting a 1 to unity. 

In the matching scheme the analytic fit form for //; on the primary ensemble 1 becomes: 

fl=c{* (l +C;P ( V] 2 ) +C{*ml + C(*m} + Cf« (m\-m m ) , (29) 

where m v = \{fh x + m y ) and we have taken advantage of the linearity of the expression to absorb 
any terms in m? into the leading coefficient. Here the superscript A{\) denotes the gauge action 
of the primary ensemble: for our choice of primary ensemble this is the Iwasaki action, labelled /. 
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The fit form describing fu on any other ensemble can then be obtained by applying equation 1251 to 
the above and replacing the a 2 coefficient with that appropriate to the particular action. 
Before continuing, it is illustrative to discuss how the a 2 coefficients for the Iwasaki+DSDR gauge 
action can be determined without having multiple lattice spacings with this action. Let us imagine 

that we have performed a global fit over the 241 and 321 ensembles as in ref. yj], and have thus 

f f f I 

determined the coefficients Cq through C3 and the Iwasaki scaling coefficient Ca . We then 

perform a fixed trajectory matching between the 321 and 32ID ensemble sets, providing us with 

Zf 2ID , Z^ 210 and F? 2ID . The fit form describing /// on the 32ID ensemble now has only one 

unknown coefficient, namely c[ K ' lD , which can be obtained by comparing any single simulated 

data point with the predicted value or by fitting over several points. In practice we would like the 

32ID data to contribute to the determination of the coefficients, thus we perform a combined fit to 

all three ensemble sets and allow c[ n,ID to be determined by minimizing the global x 2 - 

Recall that our choice of scaling trajectory defines the pion, kaon and Omega baryon masses to 

have no lattice spacing dependence up to terms &{ma 2 ) arising from the match-point dependence 

of Z/, Z/ 7 and R a . These terms are neglected by our power counting, hence the fit forms for these 

quantities contain no discretization terms. For example, the form for the Omega baryon with the 

analytic ansatz is: 

m h hh = C™ n + C" ! "m, + C™" (m y - m m ) + C™ Q (m h - m h0 ) . (30) 

The remaining analytic and ChPT fit forms can be found in section V-B of ref. yj]. Note that as we 
now measure the strange-quark dependence in the global fit rather than linearly interpolating to 
the physical strange mass prior to fitting, we include additional parameters for the heavy valence- 
quark dependence (where appropriate) and the heavy sea-quark dependence, in this order. For the 
analytic fit forms these coefficients are labelled following the existing sequence, for example the 
heavy valence and sea quark dependences of are C™ n and C™ a respectively. For the ChPT fit 
forms we label the parameters cq™ and CQ >mh for the valence and sea dependence of the quantity 
Q respectively. 

We perform our fits with the strange-quark mass expansion point m/,o set initially to the un- 
reweighted strange sea-quark mass on the 32 3 DWF+I ensemble set. This is then corrected to 
the physical strange quark mass a posteriori; with our power counting this requires only a redef- 
inition of the leading order coefficient (e.g. Cq -). For the ChPT forms we must also adjust the 
LECs in order to absorb the effect of adjusting the chiral scale to the conventional 1 GeV once 
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the lattice scale has been determined. 

Once the fits have been performed, we determine the physical up/down and strange quark masses 
(normalized to the units of the 321 primary ensemble) by numerically adjusting the quark masses 
in our fit functions such that m^/m^ and m^r/m^ match their physical values in the continuum 
limit. Here, as in ref. [1], we use m K = 135 MeV, m K = 495.7 MeV and m fi = 1672.45 MeV. The 
primary lattice spacing can then be extracted by dividing the predicted continuum value for ma 
in lattice units by its physical value. Using these results and the values of R a , Zj and Z/, found by 
fitting the data, the lattice spacings and physical quark masses for the other ensemble sets can be 
determined; we discuss this in more detail in section IVTl 

V. FIT RESULTS AND SYSTEMATIC ERROR DETERMINATION 

Following the 2010 analysis strategy, we split the chiral/continuum fits into three parts. In the first 
part, to which this section is dedicated, we performed simultaneous fits to m K , thk, mn, f K and 
over the three ensemble sets, from which we determined the physical quark masses (in matching 
scheme normalization), the lattice spacings and the scaling parameters, along with predictions for 
the physical pseudoscalar decay constants. The second set of fits were performed to Bk and the 
third to the Sommer scales ro and n; these are documented in sections IVTTl and IVrTll respectively. 
We also separate out the discussion of the determination of the physical quark masses in the ms- 
scheme into section IVTl 

A. Fit Results 

In the 2010 analysis we did not attempt to correct for finite- volume effects in our analytic fits, 
as the magnitude of the change was small with respect to the systematic error arising from the 
chiral extrapolation. However on the 32ID ensemble set we have data reaching down almost to 
the physical point, hence we might expect that the chiral systematic error will be reduced and that 
the finite-volume error may begin to dominate (as we discuss below, this does indeed seem to be 
the case). As a result, in anticipation of our later discussion, we perform our analytic fits to data 
corrected using ChPT to the infinite-volume limit. Although we do not have multiple volumes 
from which to measure the size of the correction directly, we expect that the finite- volume terms in 
NLO chiral perturbation theory will provide a somewhat reliable estimate now that we are so deep 
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in the chiral regime; we therefore estimate the finite-volume correction for each simulated data 
point as the fractional difference between the ChPTFV fit value for that point with and without 
the finite-volume terms applied. The analytic fits presented in this section are all performed to 
the finite-volume corrected data. Note that despite the near-physical pion masses on our 32ID 
ensembles, the smallest m n L is roughly 3.3, which is in fact larger than the value of 3.1 obtained 
for the lightest pion on the 321 ensembles, due to the greater physical volume of the 32ID lattice. 
As a result we do not expect our new data involving lighter quark masses to further enhance the 
finite- volume errors. 

In order to prevent accidental correlations between independent data from influencing the fit, while 
retaining the correlations between data measured on the same ensemble, we make use of the su- 
perjackknife technique to propagate the errors through our fits. A superjackknife distribution for 
a measurement is essentially a collection of independent jackknife distributions, each containing 
the fluctuations from a particular ensemble. As for the standard jackknife, any procedure, such as 
a fit or binary operation, is performed sequentially to each jackknife sample in all distributions. 
The total error on the superjackknife is obtained by evaluating the errors on each of its component 
jackknife distributions and adding these in quadrature. This technique was also used for our 2010 
analysis. 

As discussed in the previous section, our use of the strange-quark mass reweighting in the chi- 
ral/continuum fits differs from the 2010 strategy. Previously, each quantity was independently 
interpolated to the physical strange-quark mass prior to fitting; after the fit the values were up- 
dated to the new mass and the fit repeated, with this process iterated until convergence. We now 
constrain the heavy sea-quark dependence of each quantity to be the same on all ensemble sets 
and include multiple reweighted data points in the fit. As the number of reweighted masses differs 
between the ensemble sets, and considering that there are likely to be strong correlations between 
the reweighted data points, we might worry that the x 2 contributions of the data on the ensemble 
sets with more reweighted masses will be incorrectly enhanced in our uncorrelated fits. In order 
to avoid this we used only four reweighted strange-quark masses on each ensemble set, spread 
uniformly across the range. 

Upon performing the fits, we discovered significant (up to 4a) tensions between the fits and the 
pion and kaon data on the 32ID ensembles at the upper end of the reweighted mass range. How- 
ever, the upper limit of this mass range (m/, = 0.052) is considerably larger than the physical 
strange quark mass of ~ 0.047, which is actually very close to the directly simulated mass of 
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nih = 0.045. From the effective number of configurations calculated in section HH we estimated 
that reweighting to the physical strange quark mass introduces a 10-15% increase in our statisti- 
cal errors. As we go further from the simulated point we expect the accuracy of the reweighting 
procedure to further decrease due to the reduced overlap of the reweighted path integral and the 
original. At = 0.052 we found that the effective number of configurations was reduced to only 
15 on the lighter ensemble (down from 180) and 24 on the heavier (down from 148). This suggests 
that the measurements at the far end of the reweighting range are dominated by only a very small 
number of configurations and are therefore unreliable. As a result, the tension we observed be- 
tween the fits and the data at the upper end of the reweighting range is likely to be an artifact of the 
reweighting procedure. With this in mind, we repeated the fits again using four reweighted masses 
this time spread only over the range beginning at the simulated strange quark mass and ending at 
the estimated physical strange quark mass. In doing so we found that the tension disappeared. 
The inclusion of the 32ID ensembles greatly enhances the mass range over which our fits are 
performed. This should reduce the systematic error on the extrapolation to the physical light quark 
mass, and also allows us to consider removing some of the heavier ensembles from the Iwasaki 
data sets which may lie near the limits of convergence of NLO chiral perturbation theory. We 
removed the 241 m/ = 0.01 ensemble and the 321 m\ = 0.008 ensemble, as well as the partially- 
quenched data points on the lighter ensembles containing quarks with these masses. In performing 
this cut, we restrict our fits to pion masses smaller than 350 MeV, where previously the upper 
bound was 420 MeV. This amounts to a ~ 30% reduction in the largest unitary light-quark mass. 
Note that this is not a straight cut on the partially-quenched pion mass as the elimination of these 
heavy ensembles also removes a number of partially-quenched pions containing these now 'heavy' 
quarks ranging down to ~ 230 MeV. 

With the cut data set we were able to obtain excellent fits using the ChPT and ChPTFV ansatze. 
For the analytic ansatz we again found excellent fits to the decay constants as well as m# and 
mjj, but for the pion mass we found a number of outlying data points on the 32ID ensembles 
that deviated from the fit by up to 4a, with the typical size of the deviation being 0(2%). These 
deviations appear to occur due to non-linearities in the light data. The fact that no corresponding 
deviations appear for the ChPT fits suggests that these non-linearities are consistent with the NLO 
chiral logarithms. However, the discrepancies are also of the size expected for NLO terms in the 
Taylor expansion that are beyond the range of our power counting, hence we cannot draw any 
strong conclusions about their nature within our modest range of masses. As the linear ansatz 
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must be locally correct around the physical point, we sought to reduce these discrepancies by 
further lowering the cut for these fits, first by eliminating the 321 m/ = 0.006 ensemble, then by 
systematically removing the data corresponding to the heaviest partially-quenched pions. The 
limit to which this bound can be pushed in our analysis is dictated by the stability of the fits and 
the necessity to retain some data on the remaining 241 ensemble such that the a 2 coefficients of the 
decay constants can be determined; the latter implies that a 240 MeV bound is the lowest that we 
can currently reach. In practice we reached a similar level of agreement between the analytic fits 
and the data as found in the ChPT fits by lowering the bound to 260 MeV. Although this removes 
a large amount of data, we found that the fit remained very stable and that the effect of the cut on 
the values and precision of the fit parameters and predictions was surprisingly small; the typical 
change was of the order of a few percent, with the only large, statistically significant change being 
a 15% increase in the valence light-quark dependence of f n . With this in mind, we chose the 260 
MeV cut for our analytic fits. The mass combinations of the data points remaining after performing 
this cut are listed in table IXll 

As a side note, we also repeated the ChPT and ChPTFV fits to the full data set, for which the 
upper bound on the pion mass is 420 MeV. We found that, even over this large range, the NLO 
SU(2) ChPT fits were able to describe all of our data with only a few points on the 32ID ensembles 
deviating by between 2 and 3a. 

In table IXIII we list the results for the inverse lattice spacings and quark masses obtained using each 
fit ansatz, alongside the associated uncorrelated ^ 2 /dof. The results are completely consistent, 
which suggests that the extrapolation to the physical quark masses is under control. A similar 
degree of consistency can be seen between the fit parameters (where applicable) given in table IXIIH 
Here, as mentioned in the previous section, we have adjusted the chiral scale A^ of the ChPT LECs 
to the conventional 1 GeV In figures [13] and [14] we overlay our simulated data for m n and f n on 
the 32ID ensembles with the ChPTFV and analytic fit curves respectively, and in figure [L5] we 
present similar plots for m# and fx overlaid with the ChPTFV fit curves. We list the individual 
predictions for f K , fx and their ratio at the simulated lattice spacings and the continuum limit in 
table IXIVI In figure [J6] we plot the chiral extrapolation of f % overlaying the data corrected to the 
continuum limit. 

The uncorrelated ^ 2 /dof are all less than unity, suggesting that the fits are behaving. In order to 
demonstrate the quality of the fits in greater detail, we present histograms of the deviation of the 
fit from the data in units of the statistical error in figure [T2l 
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Ensemble set 


mi 


{ m x} 


{my} 




0.008 






321 


0.006 








0.004 


0.002 


0.002, 0.004 


241 


0.01 
0.005 


0.001 


0.001 


32ID 


0.0042 
0.001 


0.0001,0.001,0.0042 
0.0001,0.001,0.0042 


0.0001,0.001,0.0042,0.008 1 

Hexcl. [0.0042,0.008]) 

0.0001,0.001,0.0042,0.008 J 



TABLE XI. Sea and valence quark masses of the data included in the analytic fit with a 260 MeV cut on the 
pion mass. The third and fourth columns give the set of partially-quenched valence quark masses; the mass 
combinations of light-light quantities {m n and f K ) are found by combining each choice of m x with each 
choice of m y from the appropriate columns, with the exception of the [m x ,m y ] = [0.0042,0.008] points on 
the 321 ensemble set, for which the partially-quenched pion masses are above the cut. For heavy-light data 
(rriK, fa) the light valence-quarks are chosen from the {m x } column, and the heavy valence-quarks from 
the full set of simulated heavy-quark values. For niQ and the Sommer scales, all data are included on those 
ensembles not marked with a dash (-). 

In the remainder of this section, we discuss how we combine the results of our fits into predictions 
for f % and fx and final results for the lattice spacings and physical quark masses (in the matching 
scheme). 
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Analytic ChPT ChPTFV 


£ 2 /dof(32IW) 


0.279(64) 0.191(55) 0.221(57) 


ami(32l) 
am s (32l) 
a' 1 (321) 


0.000320(42) 0.000307(34) 0.000308(35) 
0.02660(98) 0.02650(85) 0.02627(89) 
2.295(40) GeV 2.302(35) GeV 2.310(37) GeV 


ami(24l) 
am s (2A\) 
a.- 1 {241) 


-0.001754(83) -0.001757(75) -0.001749(78) 
0.0337(18) 0.0338(13) 0.0336(13) 
1.743(43) GeV 1.743(30) GeV 1.747(31) GeV 


am t (32TD) 
am s (32lD) 
a-\32lD) 


-0.000090(34) -0.000096(21) -0.000090(22) 
0.04667(76) 0.04674(60) 0.04671(61) 
1.372(10) GeV 1.371(8) GeV 1.371(8) GeV 



TABLE XII. The # 2 /dof, unrenormalised physical quark masses in bare lattice units (without m lcs included) 
and the values of the inverse lattice spacing aT x obtained by fitting to data with m n < 350 . 



Parameter 



ChPT 



ChPTFV 



Parameter 



Analytic 



7' 

4 
K 



R'j, 



0.983(14) 
0.929(15) 
0.9730(94) 
0.939(13) 
0.7571(65) 
0.5955(72) 



0.981(14) 
0.930(15) 
0.9719(95) 
0.935(13) 
0.7562(66) 
0.5934(76) 



0.992(21) 
0.936(16) 
0.976(14) 
0.940(14) 
0.7595(90) 
0.5976(74) 



B 



(2) 
(2) 



f 

r(2) 
^5 

r( 2 ) 



4.174(83) 
0.000616(22) 
-0.000131(69) 
-2.1(2.5) 
0.1167(31) 
-0.021(70) 
0.040(45) 
0.000560(51) 
-0.00014(13) 



4.148(86) 
0.000610(23) 
-0.000159(72) 
-2.5(2.5) 
0.1196(31) 
-0.031(68) 
0.014(43) 
0.000524(51) 
-0.00020(14) 



^0 

rf*>' 



c 



C 



0.00043(23) 
7.70(16) 
0.173(40) 
-0.041(26) 
0.1221(30) 
-0.064(88) 
0.030(47) 
1.054(32) 
0.88(16) 
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Parameter 


ChPT 


ChPTFV 


Parameter 


Analytic 


C f x , mi, 


0.422(90) 


0.484(89) 


c{ n 


0.120(49) 




0.2365(77) 


0.2364(80) 





0.2364(88) 


h 


0.01907(77) 


0.02028(76) 


C'" K 


3.637(99) 




0.00220(75) 


0.00233(80) 




0.47(20) 




3.811(61) 


3.828(64) 


c m K 


3.802(71) 




0.033(43) 


0.031(43) 


C' nK 


0.001(64) 


/(AT) 


0.1466(36) 


0.1484(37) 


r<fK 



0.1500(35) 


^ Wxl 
/ l ' 


-0.034(57) 


-0.040(57) 


r h,i 
'-a 


-0.075(69) 


„/o 


0.020(38) 


0.008(38) 


'-a 


0.013(38) 




0.00622(22) 


0.00601(22) 


c{ K 


0.349(44) 


h 


-0.0034(19) 


-0.0032(20) 


c[ K 


0.76(19) 




0.2917(42) 


0.2923(42) 


c{ K 


0.2967(64) 


^ ' fK-, m h 


0.118(40) 


0.118(40) 


c{ K 


0.144(57) 




1.6659(100) 


1.666(11) 





1.6657(99) 




2.9(1.2) 


3.1(1.2) 




3.0(1.8) 


c m fl ,m v 


5.439(58) 


5.462(63) 




5.441(65) 


c m a ,m h 


0.74(29) 


0.87(31) 




0.35(39) 



TABLE XIII: The fit parameters of each of our chiral ansatze obtained by fitting to data with m n < 350 
MeV. The parameters are given in GeV" for the appropriate power of n, and with the heavy quark mass 
expansion point adjusted to the physical strange quark mass. We have ordered the table such that the 
equivalent parameters of the ChPT and analytic fits lie on the same line. The coefficients of the chiral 
logarithms have also been adjusted so that they are defined at the conventional chiral scale A x = 1 GeV. 



B. Combining Results and Estimating Systematic Errors 

The method of combining the results obtained using our three chiral ansatze into a final prediction 
was discussed at length in ref. yj]. The main issues were firstly deciding which result or combi- 
nation of results to use for the central value and secondly deciding how to estimate the systematic 
errors arising from finite-volume corrections and the extrapolation to the physical quark masses. 
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f32l 
J 7t 

pA\ 
J 71 

fill 
J 71 

J 71 



Analytic ChPT 



ChPTFV 



0.1249(21) 0.1242(22) 0.1264(22) 

0.1238(28) 0.1238(25) 0.1258(26) 

0.1284(18) 0.1273(18) 0.1280(18) 

0.1264(28) 0.1247(27) 0.1271(27) 



fJUW 
JK 

f2<UW 
JK 

fillD 

JK 

/■continuum 
JK 



Analytic ChPT 



ChPTFV 



0.1502(23) 0.1499(23) 0.1512(24) 

0.1485(30) 0.1491(26) 0.1503(27) 

0.1536(21) 0.1526(21) 0.1531(21) 

0.1525(28) 0.1509(29) 0.1524(30) 





Analytic ChPT 


ChPTFV 


(fx//*) 32 ™ 


1.202(12) 1.207(9) 


1.197(9) 




1.199(18) 1.205(11) 


1.195(11) 


[hlfnT D 


1.196(4) 1.199(4) 


1.196(4) 


( f If \ continuum 


1.206(14) 1.211(12) 


1.199(12) 



TABLE XIV Predictions for f % (top-left) and (top-right) in GeV as well as their ratio (bottom) for each 
global fit ansatz at each simulated lattice spacing and in the continuum limit obtained by fitting to data with 
m K < 350 . 



The discussion was focussed on the predicted decay constants as they are known to high preci- 
sion. We observed that our predicted value for f n from the ChPTFV fit was around 10% too low, 
and 5% in the analytic case, with smaller discrepancies visible in the kaon decay constants. We 
concluded that these are of the size expected for NNLO terms in the chiral expansion, as obtained 
by squaring the difference between our data and / - the leading order term in the ChPT chiral ex- 
pansion. Noting that both the analytic fits and ChPTFV fits appeared to describe our data equally 
well, we decided to average the two results and take their full difference as our estimate for the 
chiral extrapolation systematic. We estimated the size of the finite-volume systematic error from 
the full difference of the ChPTFV and ChPT results. 

Now that we have data ranging down almost to the physical point, we are able to revisit the issue 
of estimating the systematic errors. We first note that the differences between the ChPTFV and 
analytic results for f n and fx are now very small, smaller in fact than the formerly sub-dominant 
finite-volume contributions estimated from the difference between the ChPT and ChPTFV results. 
By comparing the above results with those obtained by fitting to all available data we observed 
that this reduction is mainly due to our removal of the data corresponding to heavier pion masses 
from the fits. 

As discussed in the previous section, we performed our analytic fits to finite-volume corrected data 
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(Y-Y & )/a Y (Y — y fit ) / cry 



FIG. 12. Histograms of the deviation of the fit from the data for each quantity on each of the three ensemble 
sets (321 top, 241 middle and 32ID bottom) with the analytic (left) and ChPTFV (right) ansatze. 
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in anticipation of the increased importance of these effects on our results. Here we investigate how 
large an effect the finite- volume corrections have on the analytic fits by repeating the latter with 
un-corrected data. The resulting fit parameters and predictions are compared to the original fits in 
table IXVl In the table we also provide the superjackknife ratios of the fit results with and without 
finite- volume corrections. We notice that in taking the ratio, many of the correlated fluctuations 
cancel, exposing underlying changes that were formerly masked by the statistical error. We ob- 
serve that in many cases the deviation of the ratio from unity is statistically significant but is only 
0(2%) or less; these changes are of the order expected for higher-order (mass-squared or a 2 m) 
effects that are beyond the range of our power counting, hence we cannot draw any conclusions 
from these results. The only quantities that change significantly are the slopes of f K , fx and m K 
with respect to the light-quark masses; this behaviour is expected as the finite-volume corrections 
will be larger in the light quark-mass regime, in which the physical length scales are greater. We 
observe a 1.7 MeV upwards shift in the continuum prediction for f K , which is consistent with the 
2.4 MeV difference between the ChPT and ChPTFV results. 

Although we now correct for the finite-volume using NLO chiral perturbation theory, we note that 
resummation techniques [24] may lead to somewhat larger estimates of the finite-volume effects. 
As we lack the ability to repeat our calculations on a larger volume, we choose to continue to 
include a conservative finite-volume systematic error in our final results, obtained, as before, from 
the full difference of the ChPTFV and ChPT results. 

In the previous section we demonstrated that the ChPTFV fit forms describe our data reliably 
over a considerably larger range of pion masses than the linear ansatz. For the final predictions 
given in the following sections we therefore take the ChPTFV results for our central values and 
use the analytic ansatz only to estimate the chiral systematic. However, we continue to find it 
surprising that a linear ansatz appears capable of describing QCD at the 1% level from the 260 
MeV pion-mass regime down to the physical point, and at the 2% level if that range is extended to 
350 MeV. 

In some cases we observed that the superjackknife errors on the differences between results ob- 
tained using the three parameterizations were larger than the differences between the central val- 
ues. In these cases we chose to be conservative and took the statistical error on the difference for 
our estimate of the systematic error. 



Quantity 


Original data 


FV corrected data 


Ratio 7? 


\R-l\/a 


£ 2 /dof 


0.219(54) 


0.274(65) 


1.254(33) 


7.672 


mi 


0.002230(58) 


0.002259(58) 


1.01293(53) 


24.379 


m h 


0.0627(12) 


0.0626(12) 


0.99857(31) 


4.603 


Zi(7AI) 


0.996(22) 


0.992(21) 


0.99619(41) 


9.377 


Zi(32ID) 


0.927(16) 


0.936(16) 


1.00932(42) 


22.131 


Z h {2AI) 


0.975(14) 


0.976(14) 


1.00073(20) 


3.580 


Z h (32ID) 


0.942(14) 


0.940(14) 


0.99876(29) 


4.366 


R a {24l) 


0.7595(91) 


0.7595(90) 


1.00004(17) 


0.218 


R a {32ID) 


0.5977(75) 


0.5976(74) 


0.99980(22) 


0.911 


cr x (327) 


2.295(40) 


2.295(40) 


1.00025(24) 


1.032 


a" 1 (247) 


1.743(43) 


1.743(43) 


1.00029(41) 


0.710 


a~ l (32ID) 


1.372(10) 


1.372(10) 


1.000048(21) 


2.282 


fn 


0.1247(27) 


0.1264(28) 


1.01359(72) 


18.879 


h 


0.1515(28) 


0.1525(28) 


1.00627(31) 


19.916 


fc/fn 


1.213(12) 


1.202(12) 


0.99143(36) 


24.076 




-0.00011(16) 


-0.00014(16) 


1.28(51) 


0.563 




3.378(30) 


3.355(30) 


0.99334(28) 


24.172 


C mic 


0.084(18) 


0.075(18) 


0.892(20) 


5.490 


3 


-0.016(11) 


-0.018(11) 


1.084(74) 


1.146 




0.0539(13) 


0.0547(13) 


1.01534(81) 


18.970 


'-a 


-0.013(17) 


-0.012(17) 


0.91(13) 


0.702 




0.0093(91) 


0.0057(89) 


0.61(38) 


1.032 


cf* 

1 


1.121(32) 


1.054(32) 


0.9404(21) 


28.102 




0.94(16) 


0.88(16) 


0.9414(88) 


6.662 




0.120(48) 


0.120(49) 


1.001(10) 


0.110 





0.06589(63) 


0.06597(62) 


1.00113(11) 


10.556 


r m K 
L l 


1.600(30) 


1.585(30) 


0.99058(37) 


25.545 


r m K 
L 2 


0.208(86) 


0.206(85) 


0.99130(46) 


18.716 


r m K 
C 3 


1.6544(97) 


1.6561(97) 


1.00101(11) 


9.208 
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Quantity 


Original data 


FV corrected data 


Ratio R 


\R-\\/a 




-0.000(28) 


0.000(28) 


0(1500) 


0.009 


d K 

u 


0.0705(14) 


0.0710(14) 


1.00633(32) 


19.911 


r fK,I 

^ a 


-0.014(13) 


-0.014(13) 


1.011(22) 


0.510 




0.0041(73) 


0.0024(72) 


0.60(76) 


0.529 


c{ K 


0.378(44) 


0.349(44) 


0.9246(84) 


9.007 


d K 

z 


0.77(20) 


0.76(19) 


0.9793(31) 


6.752 


c{ K 


0.2965(65) 


0.2967(64) 


1.00060(22) 


2.749 


c{ K 


0.144(57) 


0.144(57) 


0.9982(45) 


0.395 





0.7992(100) 


0.7994(99) 


1.00023(11) 


2.056 


L l 


3.0(1.8) 


3.0(1.8) 


0.9937(29) 


2.158 


L 2 


5.436(65) 


5.441(65) 


1.00093(22) 


4.281 


l_ 3 


0.35(39) 


0.35(39) 


1.013(29) 


0.454 



TABLE XV: A comparison of the results of analytic fits to the simulated data and the data corrected to the 
infinite volume using the ChPTFV fit forms. The quantity in the fourth column is the jackknife ratio of the 
results, R, and the quantity in the fifth column is the statistical significance of the deviation of this ratio from 
unity. 



C. Global Fit Predictions 

Applying the procedure detailed above, we present our predictions for the pion and kaon decay 
constants: 

/* = 127. 1(2.7) (0.9) (2.5) MeV, (31) 
f K = 152.4(3.0)(0.7)(1.5)MeV, (32) 
/*//* = 1.1991(116)(69)(116). (33) 

Here the errors are statistical, chiral and finite-volume respectively. Note that by restricting the 
ChPTFV fit to m K < 350 MeV rather than m K < 420 MeV used in the 2010 analysis (a 30% 
cut in the light quark mass), we obtain a value for f K that is now consistent with the known 
physical value, justifying our assertion that the previously observed deviation was mainly due to 
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the influence of higher order terms in the chiral expansion. 
For the inverse lattice spacings we obtain: 



a" 1 (321) = 2.310(37) (17) (9) GeV, 



(34) 



a-\2A\) = 1.747(31)(24)(4) GeV, 



(35) 



a _1 (32ID) = 1.3709(84) (56) (3) GeV. 



(36) 



For comparison, in the 2010 analysis we obtained a" 1 (321) = 2.282(28) ( 1)( 1) GeV and a -1 (241) = 
1. 730(25) (1)(0) GeV by fitting only to the Iwasaki data. These results are statistically consistent, 
although we find a considerable enhancement in the systematic errors. Upon detailed investigation 
we determined that these differences arise almost entirely because the scaling factors Z/, Z/ 2 and 
R a are now allowed to vary between the fits (generic scaling), as opposed being fixed to the values 
obtained at some unphysical mass point (fixed trajectory) as in the 2010 analysis: In the fixed 
trajectory case the prediction for the physical Omega baryon mass, which we use to set the overall 
scale, can vary only through the values of the physical light and strange quark masses, whereas in 
the generic scaling case the scaling parameters are those that contribute to the minimisation of the 
global x 2 , an d can thus introduce larger variations in the predicted Omega mass. This does not, 
however, suggest that generic scaling is worse than the fixed trajectory approach, as the shifts in 
the scaling parameters between the three ansatze in the former approach would simply be absorbed 
elsewhere in the latter, increasing the systematic error on some other quantities. 
Using the NLO SU(2) ChPT fits we can obtain values for the effective couplings Z3 and I4. For the 
ChPTFV and ChPT fits on their own, we find: 



As before we take the ChPTFV result for our central value. Although we cannot obtain a chiral 
extrapolation error without a corresponding analytic fit result, we can continue to estimate a finite- 
volume error from the difference between the two ChPT results. Therefore, our final values for the 
effective couplings are as follows: 



where the errors are statistical and finite- volume respectively. In the 2010 analysis (applying 
the same procedure to obtain the finite-volume error), we found I3 = 2.57(18)(25) and I4 = 



h = 2.91(23), l A = 3.99(16) (ChPTFV) 
h = 2.98(22), l A = 3.90(16) (ChPT) . 



(37) 



f 3 = 2.91(23)(7), / 4 = 3.99(16)(9), 



(38) 
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3.83(9) (7). Comparing to our fits without the reduced pion mass cuts, we determined that the 
inflation of the statistical error and the rises in the central values over the 2010 analysis results 
derive mostly from the lowering of the cut from 420 MeV to 350 MeV. However the values for It, 
and I4 agree more closely in our current analysis even without the reduced cut, suggesting that the 
inclusion of the 32ID ensembles has some stabilizing influence upon the fit. For comparison, the 
FLAG working group obtained Il25ll an estimate of I3 = 3.2(8), which was chosen to cover a large 
number of independent lattice results for this quantity, among which there are some discrepancies 
between the values. Our result is entirely consistent with this estimate. For I4, the inconsistencies 
between the results were considered too large to make a meaningful estimate. For both of these 



26, 



quantities, recent results include 2+ If determinations by the MILC collaboration 
2010 analysis paper [1], and a 2+1+lf determination by the ETM collaboration [|28J]. 
Finally we give our predictions for the physical quark masses on the primary ensemble set: 



27D and our 



m ud {32l) = 2.243 (46) (24) (10) MeV, m,(32I) = 62.2(1.1) (0.5) (0.3) MeV 



(39) 



In the 2010 analysis we obtained m ud (32l) = 2.355(81)(79)(42) and m,(32I) = 63.7(9)(1)(1). 
These numbers are again consistent, although here it appears that the enhanced control over the 
chiral extrapolation afforded by the 32ID ensembles has decreased the statistical error on the aver- 
age up/down quark mass in spite of our exclusion of a large number of data points. We also observe 
a vastly improved chiral extrapolation systematic and a substantially reduced finite-volume error 
on this quantity. In the next section we discuss how these masses are renormalized into the ms 
scheme. 



VI. PHYSICAL RESULTS FOR THE LIGHT- AND HEAVY-QUARK MASSES 
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FIG. 13. Global fits obtained using NLO SU(2) chiral perturbation theory with finite- volume corrections 
for the pion mass (top) and f n (bottom) on the 32ID ensembles. Here the left-hand plot of each pair show 
the data at the simulated strange-quark mass and the corresponding fit curves on the m/ = 0.001 ensemble, 
and the right-hand plots those on the m\ = 0.0042 ensemble. The plots of the pion mass have m 2 n / (m x + fh y ) 
on the ordinate axis, a quantity used traditionally to emphasize the chiral curvature of the data. 
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FIG. 14. Global fits using the analytic ansatz with finite- volume corrected data for the pion mass (top) 
and f n (bottom) on the 32ID ensembles. Here the left-hand plot of each pair show the data at the simulated 
strange-quark mass and the corresponding fit curves on the mi = 0.001 ensemble, and the right-hand plots 
those on the mi = 0.0042 ensemble. The plots of the pion mass have m^/(m x + m y ) on the ordinate axis, 
a quantity used traditionally to emphasize the chiral curvature of the data. The circular points are those 
included in the fit, and the diamond points those excluded by the cut on data with m n > 260 MeV. 
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FIG. 15. Global fits obtained using NLO SU(2) chiral perturbation theory with finite- volume corrections 
for the square of the kaon mass (left) and fx (right) on the 32ID ensembles. 
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FIG 16. The chiral extrapolation of the pion decay constant using the analytic and ChPTFV ansatze. 
Overlaying these curves we have plotted the unitary data extrapolated to the continuum limit using the a 2 
dependence of our fit forms. The lighter-shaded points were corrected using the analytic fit form, and the 
darker points by the ChPTFV form. Here the circular points are those included in the fit, and the diamond 
points are those excluded by the cuts at 350 MeV (ChPTFV) and 260 MeV (analytic). The upper and lower 
square points show the continuum predictions obtained using the ChPTFV and analytic ansatze respectively. 
Note that the analytic fit does not include any unitary data points on the 321 and 241 ensembles as they lie 
above the pion mass cut (cf. table IXII) . 
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In the previous section we determined the physical quark masses in lattice units in the matching 
scheme defined in section[V] In this section we discuss how we convert these into the conventional 
ms- scheme. 



A. Non-perturbative Renormalization for the Quark Masses 

We cannot simulate with a non-integer number of dimensions, hence we must match our lattice 
results to perturbation theory in order to quote a result in the Ms-scheme. Rather than matching 
using lattice perturbation theory, which is often poorly convergent, we obtain the renormalization 
coefficients Z^ s non-perturbatively at each lattice spacing via several intermediate renormaliza- 
tion schemes - the so-called RI/SMOM schemes - that are variants of the Rome-Southampton 
RI/MOM scheme. In these schemes the renormalization coefficients are calculated by fixing the 
values of appropriate amputated vertex functions, constructed using quark propagators on Landau- 
gauge fixed configurations, at a renormalization scale defined by the quark momenta. These 
schemes are defined without reference to a particular regularization, hence they can easily be 
formulated in continuum perturbation theory with dimensional regularization, and the matching 
coefficients between them and the ms scheme can be determined without reference to the lattice 
regularization. The matching is performed at a sufficiently high energy scale to be within the 
perturbative regime. 

We have shown [1] that renormalizing at 3 GeV rather than the conventional 2 GeV results in 
a significant improvement in the contribution to the systematic error from the truncation of the 
perturbative series. The quark masses in ref. [1] were calculated at the 2 GeV scale; in this analysis 



we update the procedure to use the higher scale, and use twisted boundary conditions to gain better 



control of the discretization effects on the off-shell amplitudes entering the renormalization [2911 



320. 



The RI/SMOM— > ms matching coefficients at one-loop [30J] and two-loops are known [|3 1L 
For the lattice calculation of the RI/SMOM renormalization coefficients, we are constrained in our 
choice of renormalization scale only by the desire to avoid large discretization and finite-volume 
effects. Therefore for a lattice of spatial extent L and lattice spacing a, we must choose a scale /i 
in the window: 

Zr 2 <|i 2 <(^/a) 2 - (40) 
However, if we wish to match to the ms scheme, this window is further constrained to the typically 
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much smaller regime in which both the discretization and non-perturbative effects are small: 

Aqcd«M 2 «(^/«) 2 - (41) 



This is known as the Rome-Southampton window |33|]. For the 321 and 241 lattices, with a~ l ~ 2.3 
GeV and 1.75 GeV respectively, our target of 3 GeV is accessible directly within this window. 
However, for the 32ID lattice, with a -1 ~ 1.37 GeV, we cannot calculate the lattice renormaliza- 
tion conditions in the perturbative regime without incurring large discretization errors. The 32ID 
renormalization factors are not needed for the analysis of the quark masses in this section (see 
below), but this is an issue for B^; we discuss this further in section [VTT1 

The need to calculate the RI/SMOM coefficients within the perturbative regime can be circum- 



vented via the use of off-shell step-scaling functions |29|, |34J] determined through a continuum 
extrapolation of the scale dependence (with a fixed lattice action) - in this limit the dependence 
on the action disappears and the scale dependence becomes universal. Similar step-scaling func- 
tions were used in our recent analysis of the K — > %% AI = 3/2 amplitudes [7]. In that analysis, 
performed only on the 32ID ensemble set, we used the following strategy: 

1 . We evaluated the Z-factors (or the matrix of Z-factors in the case of operator-mixing) at 
a low energy scale /io on the 32ID lattice and computed the relevant renormalized matrix 
elements. The scale /io was chosen within the region given in equation |40l in which the 
finite-volume and discretization effects are small. In practice we chose /io ~ 1 . 1 GeV. 

2. We computed the scale evolution between /io ~ 1 . 1 GeV and /i = 3 GeV of these operators 
on the finer Iwasaki (IW) lattices, upon which the high scale lies within the usual Rome- 
Southampton window. At finite lattice spacing arw, if Z s (/i,arw) is the renormalization 
factor of the operator under consideration in a (lattice) scheme S, the corresponding scale 
evolution is given by 

£ S (Ai 5 Aio,<3rw) = Z s (/i,aiw) (Z s (jUo,arw)) _ • (42) 

The result was extrapolated to the continuum limit, giving the universal running in this 
energy range for this given scheme S: 

a s (/i,/i ) = lim E s (/i 7i Uo,aiw). (43) 
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3. We multiplied the Z-factors obtained in step 1 at the scale /lo by the continuum non- 
perturbative running obtained in step 2 to obtain the desired Z-factors at 3 GeV. We then 



converted these to the ms scheme using one-loop perturbation theory 13511 



Further details of the renormalization strategy used in the aforementioned analysis can be found 



inref. 113611 . 

It is conceptually cleaner to divide our determination of the Ms-scheme quark masses in a similar 
way to the above, separating the calculation of the non-perturbative renormalization coefficients 
and their subsequent continuum extrapolation from the perturbative matching stage. We therefore 
first calculate the RI/SMOM coefficients at a low energy scale /lo and then calculate the step- 
scaling functions from this scale to 3 GeV. As discussed in section Ivnl the choice /lo = 1.4 GeV 
is optimal for the Bk analysis - we use this scale in the quark-mass analysis for consistency. 
Providing the jackknife/bootstrap errors are propagated correctly, the value of Z^(3 GeV) obtained 
after applying the step-scaling function to the 1.4 GeV result will be exactly the same as if we had 
performed the continuum extrapolation directly at 3 GeV, due to the fact that the step- scaling 
functions are calculated using the same data. 

1. Determination of the Lattice Renormalization Coefficients 

Before presenting the results of our analysis, we summarize our measurement strategy, highlight- 
ing several important improvements over the original PJ/MOIV \ methods. 

The original RI/MOM scheme, defined in ref. 113 311 . was shown Il37ll to suffer from greatly enhanced 
chiral symmetry breaking errors. These were found to occur due to the use of so-called exceptional 
kinematics, for which the vertex has channels along which the momentum transfer is zero; these 
allow quark and gluon loops with momenta below the spontaneous chiral symmetry breaking 
scale to exist even when the external momenta are moderately hard. The persistence of non- 
perturbative effects at high energy gives rise to large uncertainties in the perturbative matching. 
In order to avoid this problem we follow the 2010 analysis procedure in using non-exceptional 



'symmetric' kinematics [370 for which no exceptional channels exist. With these kinematics the 
non-perturbative effects fall off much faster as the virtuality is increased. 

The quark mass renormalization coefficient Z m , which is taken in product with the bare quark 
mass to obtain the renormalized quantity, is determined from the flavor non-singlet scalar and 
pseudoscalar vertex renormalization coefficients, Z$ and Zp respectively, via the relation Z m = 
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l/Zs = 1 /Zp. The equivalence of Z5 and Zp is not exact if the chiral symmetry is broken; this 
occurs due to the low-energy spontaneous chiral symmetry breaking of QCD and to a much lesser 
degree from finite-L 5 effects. With non-exceptional kinematics, the former vanishes as l/p [|38|l . 
and is therefore small at the 3 GeV scale at which we perturbatively convert to the Ms-scheme. In 
ref. [38] and during the present analysis we found that the effect of the difference between Z5 and 
Zp on our final ms scheme quark masses was considerably smaller than the error associated with 
the truncation of the perturbative series; as a result we do not need to include a systematic error 
for this effect. For the central values we arbitrarily chose to take the average of the scalar and 
pseudoscalar renormalization factors to determine Z m , as was performed in ref. 113 811 . 
The scalar and pseudoscalar vertex functions 115 and Hp were constructed at all sink locations 
of two quark propagators with momenta p\ and p%. The symmetric kinematics require that the 
momenta are chosen such that p\ = p\ = (p\ — P2) 2 = q 2 = — jU 2 for a renormalization scale 
/l. As before we used volume momentum source propagators as these have been shown yfl to 
significantly reduce the statistical error on the NPR coefficients. 

The renormalization conditions for the scalar and pseudoscalar vertex functions, applied at the 

Zs Zp 

scale /l in the three-flavor chiral limit, are: — As = 1 and — A5 = 1, where 

^q ^q 



A, = |ltr[n,./], A^ltr 



Up-y 5 



(44) 



Here / is the identity matrix and Z q is the wave-function renormalization factor. Z m in the non- 
exceptional schemes is thus calculated as 

Z m Z q = l -(A s + A P ). (45) 

The wave-function renormalization factor is determined from the renormalization condition on the 

Zy 

vector current: — Ay = 1, where 

Zq 

Ay = ^tr[lI v r M ] (46) 

for the vector bilinear vertex Ily^. With symmetric kinematics, the momentum transfer q 2 is 
non-zero, hence we have two choices for the projection matrix T^, namely 7^/4 and $q^/q 2 ; 
these define two different renormalization schemes which we label RI/SMOM^ and RI/SMOM 
respectively. Here we have used q^ = sin(^) following ref. yj]. In the remainder of this work we 
refer to the two schemes collectively as the 'SMOM schemes'. 
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In the above, the vector renormalization coefficient Zy is identical to the factor relating the four- 
dimensional vector current to the corresponding Symanzik current. In section [TIT] we discussed 
how this quantity can be calculated independently using the ratio of the local four-dimensional and 
conserved five-dimensional vector currents. (The values for this quantity on the Iwasaki ensemble 
sets were determined in ref. [1].) As Zy is known, we can combine its measurement with the ratio 

obtained from the vector- vertex renormalization condition, in order to determine Z q . 
In principle a separate measurement of Z q could be obtained using the axial-vector vertex. As was 
the case for the scalar and pseudoscalar vertex functions, this measurement can differ from that 
calculated via the vector vertex due to the residual effects of the low-energy spontaneous chiral 
symmetry breaking and small finite-L 5 effects. However, in refs. yj] and yO we found that the 
effect of the difference between the vector and axial-vector vertex functions on our final result is 
again negligable compared to the perturbative truncation error. 

As mentioned above, the renormalization conditions are applied in the three-flavor chiral limit. In 
practice we generate data on each ensemble with quark masses set equal to the dynamical light- 
quark mass; the chiral extrapolation is then performed using a linear fit over the unitary light-quark 
mass-dependence. The vertex functions are flavor-independent, hence this extrapolation also takes 
the valence strange-quark, but not the sea strange-quark, to the chiral limit. As we have only a 
single simulated dynamical strange-quark mass and reweight over only a short range, we cannot 
reliably take this final mass to the chiral limit. In refs. yj] and we estimated the effect of 
not taking the strange sea-quark to zero using the slope of the unitary light-quark extrapolation, 
reduced by a factor of two to obtain the contribution of a single flavor. For the RI/MOM scheme, 
the two-flavor mass -dependence was found to be significant, resulting in a large systematic effect 
comparable in size to the truncation systematic. However, for the RI/SMOM schemes we found 
a very benign mass-dependence that was statistically indistinguishable from zero. Note that this 
estimate is highly conservative as the slope is likely dominated by the valence mass dependence; 
this suggests that we can ignore this systematic effect in our present analysis, for which we use 
only the non-exceptional schemes. 

In ref. yj] we calculated the renormalization factors over a range of momentum scales. The scales 
at which we could perform our lattice measurements were limited by the need to form a symmetric 
momentum configuration with spatial momentum components that are discretized in units oiln/L 
by the periodic boundary conditions. The resulting momentum configurations were typically dis- 
tinct under the hypercubic group, hence the measurements were susceptible to lattice artifacts that 
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vary under 0(4) rotations. These induced a scatter in the data, breaking the expected smooth scale 
dependence; as a result we were forced to artificially inflate our errors by a factor of ^/x 2 /dof, 
taken from a straight-line fit to the data. In ref. I129n we showed that the scatter can be eliminated 
entirely using twisted boundary conditions to induce quark momenta with a fixed direction; the 
remaining lattice artifacts can be removed by a continuum extrapolation. This approach was used 
for the renormalization of Bk in the second of the 2010 analysis papers yfl. In the current analysis 
we also adopt this technique for the quark mass renormalization. 



2. Perturbative Matching to the MS Scheme 



The conversion factors (^/ SM0M- ^ MS between the RI/SMOM and ms schemes were first computed 



at one loop in ref. 



known from refs. [31 



30 



and the two-loop corrections for both RI/SMOM and RI/SMOM r are 



32TI . Regarding our notation, we write the running of the renormalized 



quark mass (in a given scheme S) between the scale /Iq and ji in the form 



m s ( A i)=m s ( A i )expf r (At) dx|^) , 
\Ja s (no) pyx) J 



(47) 



where, following H39fl . we use a s — (a s /n). We expand the anomalous dimension y% and the 
/3 -function (dropping the superscript S for clarity) 



(>1 



P(a s ) 



(48) 
(49) 



i>0 



where we have made explicit the fact that is scheme-independent (we do not discuss here the 
scheme dependence of a s , which cancels in equationl47l). We can then express the result of eqn.l47 
with the help of 

(50) 



where 



Ufa) P(x) J c s (/io) 



c s ( M ) = a,Qi) ^ ( 1 + ( i ^ ) as {p) + 



(51) 



Still following [13911 . we then define the renormalization-group-invariant (RGI) mass m by 



_ r (0) 

m= lim m s (n)a s (n) ?o . 



(52) 
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Using equation [47] and [ST] this gives 



m s (Mo) 
c s (Mo) 



™ = TsT^T v A*o ■ (53) 



In particular, since m is renormalization group invariant, we can use the ratio of c's to change 
scheme: for example, the conversion factor between a scheme Si and a scheme S2 at the scale ji 
is given by 

C s ^(a) = ° S2 ^ (54) 
The RGI mass is obtained from m (fl) using equation [53] which implies that 

^°V) = ^. (55) 



With some simple linear algebra we are able to convert the numerical results of ref. Q30Q to our 
conventions and evaluate eqn 



Finally, to obtain a s at 3 GeV in the three-flavor theory, we used the four-loop running of refs. [40] 



and [41] and took a s (mz) = 0. 1 184 I142h as an initial condition. We ran this quantity down to the 
charm mass, changing the number of flavors when crossing each threshold, obtaining a s (3 GeV) = 
0.2454. 

Putting everything together, we found 

c Ri/SMoi (3 GeV) = 3 - 1052 C 1 " °- 0825 " 00066 + ^)) = 2-8283 , (56) 
3. 1052 (1-0. 1086- 0.0147 + @(a 3 s )) = 2.7223, (57) 



c Ri/SMOM rp (3 Gey) 



= 3. 1052 (1 -0.0699- 0.0035 + ^(a^)) = 2.8773. (58) 



c MS (3GeV) 

Combining these we obtained for the conversion factors: 

ri/smom^ms = c jM 3GeV )) = 0.9830. 

ri/smom^ms c m (a s (3 GeV)) 
C m (3 GeV) = pt /smom = °- 9462 » ( 6 °) 

c RI/SMOM yp p QeV) 

which are correct to order a?. 



(59) 



With the four- loop anomalous dimension of ref. Q39Q, we obtain 



C S ^ RGI (3 GeV) =3. 1052 (1-0.0699 -0.0035 -0.0001 + 0(a*)) =2.8769. (61) 
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For completeness we apply the same procedure at a renormalization scale of /i = 2 GeV (using 
0,(2 GeV) = 0.2960): 

2.5452 + 0(a 3 s ), (62) 



c ri/smom( 2 GeV) 
1 

c RI/SMOM 7M(2GeV) 
1 



2.4218 + 0(al), (63) 



2.6017 + 0(a 3 s ), (64) 



c MS (2 GeV) 
which are again quoted to 0(a 3 s ). Thus we find 

c RI/SMOM^MS (2 GeV) = Q 97g3 ? (65) 
RI/SMOM v -^MS 

C m 7M (2 GeV) =0.9309, (66) 

c ms->rgi^ 2 GeV ) = 2.6012 , (67) 

where, as before, the SMOM to ms conversion factors are correct up to terms 0(a 3 ), whereas 
the RGI conversion factor is true up to terms &(d%) by virtue of using the four-loop anomalous 



dimension. As expected, these numbers are in very nice agreement with the ones given in ref. [31]. 
We close this paragraph with a remark about the definition of the RGI quantities: Our convention 
is such that, at the first order of perturbation theory, the conversion to the RGI quark mass is given 
by 

cr RG V) = (aM/x)- 4/{n - 2llf/3) ■ (68) 
This convention differs from the one used for , where a s is not divided by n: 

C^ RGI (A0 = (a,(Ai)n 2/(11 - 2 "/ /3) ■ (69) 

(Here the difference between the anomalous dimensions of the quark mass and Bk accounts for the 
factor of two in the exponent between the two expressions.) Although the difference in conventions 
is rather unfortunate, we adopt them in order to match those commonly used in the literature. 



3. Calculation ofZ m in the RI/SMOM Schemes 

We calculated the RI/SMOM and RI/SMOM yM bilinear vertex functions on each ensemble of the 
Iwasaki lattices, using quark propagators with the (twisted) momenta given in the first two blocks 
of table IXVIl (as explained below, the 32ID renormalization factors are not needed). These were 
then linearly extrapolated to the two-flavor chiral limit. We plot the scale dependence of the 



60 



241 


Pi Pi 


o 




(-2,0,2,0) (0,2,2,0) 


{-0.45136,0.732} 




(-3,0,3,0) (0,3,3,0) 


^n: ii = {-2,1..., 12} 




(-4,0,4,0) (0,4,4,0) 


3 
2 


321 


Pi Pi 


e 




(-2,0,2,0) (0,2,2,0) 


{-0.413,0.783} 




(-3,0,3,0) (0,3,3,0) 


l 

4 




(-4,0,4,0) (0,4,4,0) 


H.|> 




(-5,0,5,0) (0,5,5,0) 


I 8 ' 8 J 


32ID 









(-3,0,3,0) (0,3,3,0) 


{0.0} 




(-4,0,4,0) (0,4,4,0) 


\n : n = {-1,0,1,2} 




(-5,0,5,0) (0,5,5,0) 


\n : n = {-l,0} 



TABLE XVI. Non-exceptional momenta and twist angles used for the evaluation of amputated twisted 
Green's functions in our NPR analyses. The momenta here are listed in (x,y,z,t) order. The integer Fourier 
mode numbers {«;} are related to the lattice momenta via ap, = The momentum added by the twist is 
determined by the twist angle d giving ap; = ( 2 ' 1 '+ 9 ^ 7r . The twists that are not multiples of | are chosen to 
match specific momenta on a larger volume lattice that will be described in a forthcoming publication. 

resulting chiral-limit renormalization factors in figure \T7\ here we clearly see the smooth scale 
dependence arising from the use of twisted boundary conditions. 

In order to obtain the values at 1 .4 GeV we performed an interpolation over several data points 
in the region surrounding the 1.4 GeV renormalization scale. Using the lattice spacings from 
section |V] we find the corresponding values of (ap) 2 to be 0.367 and 0.642 on the 321 and 241 
lattices respectively. In this region of figure [FT] we see a non-linear scale dependence arising 
from the (supressed) poles and the renormalization group running, hence we cannot perform our 
interpolation using a simple linear function. Upon experimenting with several different non-linear 
forms, we found that the following parameterization: 

Z m [(ap) 2 ]=C + C l /(ap) 2 + C 2 (ap)\ (70) 

fit the data well and was stable when the number of points was increased. We present the results 
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FIG. 17. Z,„ in the two SMOM schemes at a range of scales on the 241 (left) and 321 (right) ensemble sets. 



4(/i = 1.4GeV) 
Z,^ =3.0 GeV) 



321 

RI/SMOM RI/SMOM 



241 

RI/SMOM RI/SMOM 



1.7763(43) 1.9558(36) 
1.4579(2) 1.5419(2) 



1.7782(62) 1.9612(52) 
1.4414(5) 1.5183(2) 

TABLE XVII. Renormalization factors in the intermediate RI/SMOM scheme S at the scale jU. Here the 
quoted error contains only the statistical contributions from the amputated vertex functions, not the fluctua- 
tions from the uncertainties on the lattice spacings and Zy. 



of interpolating to ji = 1.4 GeV in table IXVlTl In order to later obtain the step-scaling factors, we 
repeated the above with a 3.0 GeV renormalization scale; these results are also given in the table. 
Note that the error quoted for the results in this table contains only the statistical contributions 
from the amputated vertex functions; the fluctuations arising from the statistical and systematic 
uncertainties on the lattice spacings and Zy are discussed below. 



4. Renormalization of the Continuum Quark Masses 

The physical quark masses determined in section|V]are quoted in the 'matching scheme', whereas 
the renormalization factors above act upon the bare physical quark masses. Therefore in order to 
obtain the quark masses in either the ms scheme or one of the Rome-Southampton schemes, we 
must first convert the matching scheme masses into bare masses using equation ITTl 
The matching scheme is a non-continuum (due to its explicit cutoff dependence), mass -independent 
scheme in which a bare quark mass in physical units that is determined at a coupling /3 is renor- 
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malized by fixing its value to that obtained on a 32 3 x 64 x 16 domain wall lattice with the Iwasaki 
gauge action at /3 = 2.25: 

m match = Z match (m f ,p)m f (p) . (71) 

As discussed in section[lV] the renormalization factor Z match (m/, j8 ) at finite lattice- spacing is only 
weakly dependent upon the mass, hence we require just two factors: one to renormalize heavy 
quarks near the physical strange quark mass, and one to renormalize the light quarks. We labelled 
these Zfj(P) and Z/(/3) respectively, and calculated their values on the 241 and 32ID lattices as part 
of our global fits in section IVl (the values on the 321 ensemble are unity by definition). 
Given the values of m^Jf h and m™ atch , we can obtain quark masses renormalized in one of our 
intermediate RI/SMOM schemes S at a given /3 using the non-perturbative renormalization factors 
calculated above via the following ratios: 

m t//P) = Kft >< 4(J8)A(J8) and mf (j8) = x Z% (P)/Z h (fi) . (72) 

These quantities still retain lattice artifacts which must be removed via a continuum extrapolation. 
Since, by definition, Z/ 7 and Z/ absorb the coupling dependence of the quark masses, we need only 
extrapolate the ratios 

Z S ml (P)=Z s m (P)/Z 1 (P) and Z s mh (f3)=Zl(f3)/Z h (P). (73) 

For this we assume a linear dependence on a 2 , neglecting the higher order effects. Note that we 
cannot include the values of Z m calculated on the 32ID lattice in this extrapolation due to this 
lattice having a different gauge action, and hence a different scale dependence, than the 241 and 
321 lattices. As a result we have not analyzed this quantity in the present analysis. 
In order to correctly propagate the statistical errors and the chiral/finite-volume errors on the var- 
ious quantities we use the superjackknife procedure as before and repeated the analysis using Z/, 
Z/j, the lattice spacings and the quark masses calculated using each of the three chiral ansatze 
separately, taking the differences between these results at the final stage to determine the system- 
atic errors in the usual way. In practice, the determination of the renormalization coefficients was 
performed using bootstrap resampling and used only the final results for the lattice spacings in 
determining the renormalization scale. In order to ensure that the systematic and statistical er- 
rors were correctly propagated we devised a procedure for generating suitable 'super-jackknife' 
distributions from these; this procedure is given in Appendix |B] 
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FIG. 18. The continuum extrapolations of Z ml (left) and Z mh (right) in the RI/SMOM at 1.4 GeV. 



We performed the continuum extrapolation of Z^{ji = 1.4 GeV) and Zj, (/x = 1.4 GeV) for each 
choice of scheme S and chiral ansatz A, obtaining the values listed in table IX Villi In the table and 
below we add a superscript 'c' to denote continuum quantities. An example of the continuum 
extrapolation is shown in figure [T8l 

The step-scaling factors a s,A (3 GeV, 1.4 GeV) were then determined via a continuum extrapola- 
tion over the Iwasaki lattices of the ratio E of renormalization coefficients at 3 GeV and 1 .4 GeV 
(cf. equation 1431. This was repeated for each scheme and chiral ansatz, giving the values also 
listed in table IXVlIll We then applied the step-scaling factors to Z c ml and Z c mh at the 1.4 GeV scale 
to obtain the corresponding values at 3 GeV; these are again listed in table IX Villi Note that there is 
a quite considerable cancellation between the statistical fluctuations on the step-scaling factors and 
the 1.4 GeV renormalization coefficients; this cancellation is necessary to reproduce the smaller 
statistical errors on the 3 GeV factors and justifies the use of superjackknife error propagation. 
(Similar results might be obtained using bootstrap resampling for all quantities, with a consis- 
tent number of bootstrap samples, although this risks accidental cancellation between ostensibly 
uncorrected fluctuations.) 



5. MS-scheme Renormalization Factors and Systematic Errors 

Applying the perturbative conversion factors to Z c ml and Z c mh at 3 GeV, we finally obtain the ms 
renormalization coefficients for the quark masses determined in section [V] We list the values in 
table lXVflll All that remains prior to obtaining the ms quark masses is to decide which intermediate 
scheme to use for the renormalization and to analyze the systematic errors. 
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Ansatz A 




Quantity 


Scheme S 


Scale(s) jJ. 


ChPTFV 


ChPT 


Analytic 


7 C 


SMOM 




1.735(36) 1.735(36) 1.752(51) 


Z c 


SMOM f 
SMOM 


1.4 GeV 


1.918(39) 
1.712(27) 


1.917(39) 1.935(55) 
1.711(27) 1.712(34) 


tnh 


SMOM f 




1.893(29) 1.890(29) 1.890(37) 


a 


SMOM 

1 

SMOM f 


4 ^3.0 GeV 


0.797(8) 
0.755(7) 


798(8) 
0.756(7) 


0.799(8) 
0.758(7) 


7 C 


SMOM 




1.383(27) 


1.385(27) 


1.401(40) 


^ml 


SMOM f 




1.449(28) 


1.450(28) 


1.466(42) 




ms (via SMOM) 




1.360(26) 


1.361(26) 


1.377(40) 


7 C 


Ms (via SMOM 7 m) 
SMOM 


3.0 GeV 


1.371(26) 
1.365(18) 


1.372(26) 
1.365(18) 


1.387(40) 
1.368(25) 




SMOM f 




1.429(18) 


1.429(18) 


1.432(26) 




ms (via SMOM) 




1.341(17) 


1.342(17) 


1.345(25) 




Ms (via SMOM 7 m) 




1.352(17) 


1.353(17) 


1.355(25) 



TABLE XVIII. The factors Z c m[ and Z L mh used to convert our matching-scheme physical quark masses into 
each intermediate NPR scheme at 1.4 and 3.0 GeV, and the step-scaling factors used to run between those 
scales. We also list the ms renormalization factors with \i = 3.0 GeV, obtained by applying the perturbative 
conversion from each of the intermediate RI/SMOM schemes. The superscript 'c' on the renormalization 
factors is used to indicate that these are continuum quantities. The right-most columns correspond to the 
three choices of chiral ansatz used to obtain the lattice spacings used for the scale-setting and continuum 
extrapolations. 
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In the 2010 analysis we decided that the most reliable ms renormalization coefficients were ob- 
tained using the SMOMyP intermediate scheme. This was based on the fact that this scheme 
showed a considerably smaller scatter from 0(A) -symmetry breaking lattice artifacts than the 
SMOM scheme. However, now that the scatter has been eliminated through the use of twisted 
boundary conditions, we base our choice of 'best' scheme on the size of the error in the matching 
of the intermediate scheme to ms, which we estimate from the size of the two-loop terms in equa- 
tion [60l We see that the SMOM-scheme conversion factors appear to converge faster than those 
in the SMOMyM -scheme, with a two-loop term roughly 75% smaller. As a result we adopt the 
SMOM scheme for our final numbers. 

We expect the main contribution to the systematic error to be associated with the truncation of 
the perturbative expansion of the ms scheme-change factors. In ref. QJLD we discussed two suitable 
methods for estimating this error: The first is to use the size of the two-loop term in the perturbative 
conversion and the second to take the full difference between the ms coefficients calculated at 3 
GeV using our two intermediate SMOM schemes. For the 2010 analysis, the most conservative 
estimate was obtained from the size of the two-loop term, however, now that we have adopted 
the RI/SMOM scheme for our final result we find that the 0.4% two-loop contribution is smaller 
than the 0.8% difference between the results obtained via the SMOM and SMOMyM intermediate 
schemes. We therefore use the latter as our estimate of the truncation error. 
In section |VI A II we detailed several additional sources of error in our renormalization procedure 
that arise from non-perturbative effects; specifically, we highlighted the effects of the low-energy 
spontaneous chiral symmetry breaking and those associated with the dynamical strange sea-quark 
mass-scale. There are also likely to be additional effects at the Aqcd scale that were not con- 
sidered. Although we concluded that the non-perturbative effects at the 3 GeV matching scale 
are negligable compared to the truncation error on our final results, it is illustrative to consider 
at what point they enter into our calculations. The RI/(S)MOM schemes are actually defined in 
the limit /i 2 Ag C£) , at which the behavior is purely perturbative. The momentum schemes that 
we actually implement on our lattice can be therefore be regarded as different schemes that take 
into account the non-perturbative behavior. We therefore consider the aforementioned errors not 
as properties of the numerical renormalization factors, but rather as additional errors on the per- 
turbative conversion to the Ms-scheme, arising from the fact that the scheme-change factors are 
calculated using a slightly different scheme than the numerical results. 

There are two final sources of systematic error on the renormalization conditions - those arising 
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from the chiral extrapolation and finite-volume errors on the lattice spacings used in the scale- 
setting and the continuum extrapolation. In the previous section, we repeated the analysis using 
the lattice spacings obtained from our global fits with the three different chiral ansatze. We can 
therefore estimate these errors using the procedure discussed in section IVBl namely estimating the 
chiral systematic error as the larger of two values, the first being the difference in central values 
between the results obtained using the ChPTFV and analytic parameterizations, and the second the 
superjackknife error on this difference. The same procedure is applied to the ChPTFV and ChPT 
results to estimate the finite-volume error. We take the central value and statistical error from the 
ChPTFV ansatz. 

The final values for the quark mass renormalization factors are: 

Z c ml (m,3 GeV) = 1.360(26) (22) (2) (11) , 
Z£ A (Ms,3GeV) =1.341(17)(15)(1)(11). 
Here the errors are due to statistical, chiral, finite- volume and truncation effects. 



B. Results for the Physical Quark Masses 

Multiplying Z m \ and Z„ ; / 7 by the physical quark masses in the matching scheme, we obtain 

m lld (Ms, 3 GeV) = 3.05(8) (6) ( 1 ) (2) MeV, m s (m, 3 GeV) = 83.5(1 .7) (0.8) (0.4) (0.7) MeV, 

(75) 

where the errors are statistical, chiral, finite- volume and from the perturbative matching. The quark 
masses obtained in our 2010 analysis were quoted in the ms scheme at 2 GeV. In order to facilitate a 
comparison between these and our new results we must therefore convert to a common scheme; for 
this we use the Renormalisation-Group Invariant (RGI) scheme, for which the conversion factors 
from ms are given in eqns. [67] and [61] for 2 and 3 GeV respectively. Applying the latter to the 
results above we find: 

m ud = 8.78(24) (17) (3) (7) MeV, m, = 240. 1(4.8) (2.4) (1.2) (2.0) MeV, (76) 
where the hat is used to label the RGI values. In the 2010 analysis we obtained 

thud = 9.34(34) (31) (16) (21) MeV, m s = 250.2(3.9) (0.5) (0.3) (5.5) MeV. (77) 

Our new result appears to be consistent with that of the 2010 analysis, but has a renormalization 
systematic error that is over a factor of two smaller by virtue of performing the matching to the 
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ms scheme at 3 GeV, rather than 2 GeV, at which the perturbation theory is more reliable. For 
the up/down quark mass we also see a substantial improvement in the chiral and finite-volume 
systematics, resulting from the lowering of the pion mass cut in the fit and the inclusion of the 
32ID data. For the strange quark mass, the 32ID data does not have the same effect because 
the Iwasaki data were already (after reweighting) at the physical mass, and the light-quark mass 
dependence of the kaon is small. The larger chiral and finite-volume systematics on this quantity 
likely arise from allowing the scaling parameter Z/ ; , and also to a lesser extent Z/, to differ between 
the fit ansatze rather than remaining fixed; this allows the larger changes in the quality of the fit 
for the other fitted quantities to influence the kaon fit. A similar effect was observed for the lattice 
spacings and was discussed in section fVCl 

For comparison with the above, the FLAG working group give m M£ /(Ms,2 GeV) =3.43(11) MeV 



MILC 127, 



and ot 5 (ms, 2 GeV) = 94(3) MeV Il25h . These values were obtained by combining results from the 



4311 and HPQCD H44I1 collaborations, as well as our 2010 analysis results. Converting 
to the RGI scheme using the conversion factor given above, these become m U( i = 8.92(29) MeV 
and m uc j = 245(8) MeV, which both agree very well with our results. 

Finally, for completeness we also calculate the ratios of the strange and up/down quark masses: 



tfl 

27.36(39)(31)(22)(0), (78) 



mud 



where the errors are again as above. 



VII. CHIRAL/CONTINUUM FITS AND PHYSICAL RESULTS FOR B K 

In this section we present our results for the neutral kaon mixing parameter B^. Continuum results 
are obtained by performing chiral/continuum fits over our three ensemble sets following the strat- 
egy outlined in section [IV] This analysis extends that in ref |2|] through the inclusion of the 32ID 
ensemble set. 

As Bk is a scheme-dependent quantity we must perform our fits to renormalized data. We de- 
termine the renormalization factors again using variants of the RI/MOM scheme with symmetric 
kinematics. We first outline this calculation, then discuss the application of our chiral fitting tech- 
niques to this quantity. Finally we present the continuum results in the ms scheme at 3 GeV. 
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A. Non-perturbative Renormalization Factors 

Unlike in the case of the quark mass renormalization, we require renormalization factors for 
on both the Iwasaki and Iwasaki+DSDR ensemble sets. In this case, the option of calculating 
our lattice renormalization factors directly at 3 GeV is not an option since we cannot simulate 
within the perturbative regime without incurring large lattice artifacts. (We remind the reader that 
perturbation theory is required to match the renormalization factors computed on the lattice to a 
continuum scheme, typically MS in which the Wilson coefficients are computed.) As discussed 
in section |Vll our analysis QVU of the AI = 3/2 K — > %% amplitudes had a similar issue, which 
was solved by computing the renormalization factors at a low energy scale, )1q = 1.1 GeV, at 
which finite- volume effects and lattice artifacts are small (i.e. satisfying eqn. l40l) . and using the 
continuum step-scaling factors to evolve this to the perturbative matching scale. For this analysis 
we adopt a similar procedure. 



1. Determining the NPR factors 

We follow ref. (2j] in calculating the renormalization factors in four different lattice schemes. First 
we consider the process 

d{pi)s(-p 2 ) d(- Pl )u(p 2 ) (79) 

with a variety of momenta satisfying the symmetric momentum configuration p\ = p\ = (p\ — 
P2) 2 = /i 2 . We write the corresponding amputated Green's function evaluated on Landau gauge- 
fixed configurations as A^ U g (the color indices i, j,... and Dirac indices a, /3, . . . correspond to 
the external states). We have to project these Green's functions onto their Dirac-color structure, 
where, as before we, define two projectors using both the y-matrices and $ (where Nc is the number 
of colors and = s'm(q^)) : 



)ij,kl 



8 ij 8 kl 



1287Y c (7Y c + r 



(80) 
(81) 



These act on A in the following way: 



m = PX Al = P ij ' kl A ij ' kl 



(82) 
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As before we can renormalise the quark field, and hence obtain Z q , in both the RI/SMOM and 
RI/SMOMy/j schemes; we therefore have four independent renormalization schemes for Zb k : 

Z ( ^ ) = (zf)) 2 [pW{A}]- 1 (83) 

where A and B can be either or Here the label (27, 1) refers to the SU (3)l x SU (3)r trans- 
formation properties of the VV +AA four-quark operator that forms the numerator of equation [TBI 
Motivated by [2J], we focus only on two schemes : the (A, B) = (7^,7^) and (^, 4) combinations. 
The renormalization factor for B& is then 

7 (A,B) 

2g* = -^. (84) 

We obtain Z^jZ\ from the renormalization conditions on the vector and axial-vector vertices: 

|* = i(A A +Av). (85) 

As discussed in section [VIl the difference between these vertices in the SMOM schemes is tiny 
and can be ignored; we used their average only such that the same procedure can be applied for 
the exceptional schemes. 



2. Perturbative Conversion Factors 



The one-loop perturbative conversion factors for converting to the ms- scheme from the SMOM 
schemes are obtained using the expressions in ref. yfl, resulting in the following: 

C { Jf = 1 - 0.45465 (H) =0.99112 and 
c {f,t) =1 + 0.21197 (H) = 1.00414, 

where 

a, (3 GeV) =0.24544. (87) 

As discussed in the following section, we do not use the SMOM(^, 7") or SMOM(7^, ^) schemes 
for our final predictions, hence we have not listed the corresponding conversion factors above. 



3. Renormalization Scales 

As the 3 GeV matching scale lies within the Rome-Southampton windows for the two Iwasaki 
lattices, we need only compute the 32ID renormalization factors at the low energy scale and sub- 
sequently use the continuum step-scaling factors to run these up to the same scale as the Iwasaki 
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coefficients. However in practice we found that the statistical errors on the step-scaling factors 
were quite large, which resulted in considerably larger errors on the 3 GeV renormalization fac- 
tors than their Iwasaki counterparts. Note that contrary to the case of the mass renormalization, no 
cancellation occurs between the statistical fluctuations on Zs K (fM)) and c(3 GeV,/^) as the data 
sets from which they were determined are entirely independent. 

The disparity in the statistical errors between the renormalization factors has the effect of weaken- 
ing the constraints that the 32ID data imposes on the simultaneous chiral/continuum fit under the 
global X 2 minimization. As a naive test of the impact of this disparity, we repeated our fits with 
the errors on the 32ID renormalization factors artificially reduced to match those on the Iwasaki 
lattices. We found that the central value of the continuum prediction for Bk shifted by an amount 
comparable to the chiral and finite-volume systematics; an effect too large to be ignored. As we 
pointed out in section |V] when discussing the number of reweighting samples to use on each lattice, 
it is important to treat each ensemble set uniformly such that the weight of each of the ensemble 
sets in the fit depends only on the statistics of the data. We therefore calculate the renormalization 
factors for all three lattices at the same scale, chosen within the regime in which the discretization 
effects are under control. The 1.1 GeV scale used in ref. [7] meets this criteria, although we found 
a noticable reduction in the statistical errors by raising this to 1 .4 GeV (actually 1 .426 GeV, the 
nearest scale at which we had a simulated point). Of course, using a larger scale increases the size 
of the discretization effects on the 32ID lattice, however, as we ultimately perform a universality- 
constrained continuum extrapolation, only the &\{ap)^\ terms and higher remain in the final result 
for Bk- Only after performing the continuum limit do we apply the step- scaling factor to evolve 
the continuum prediction to 3 GeV, at which the matching to ms is performed. 



4. Results 

Following the above strategy we calculated Zg K at /lo = 1.426 GeV on each of the three ensemble 
sets. In addition, we re-calculated the Iwasaki renormalization factors at 3 GeV such that we could 
obtain the continuum step-scaling functions. The quark momenta used in these measurements are 
listed in table IXVIl and we present the values at both renormalization scales in table IXIXl We 
used the central values of the lattice spacings given in section IV CI to set the physical scales in 
these determinations. 

In order to correctly propagate errors on the lattice spacings, we formed superjackknife distribu- 
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Quantity Projector P 


ocaie jU 


Value 

321 241 32ID 




44) 




1.0608(12) 1.0320(11) 0.9992(11) 








0.9788(9) 0.9527(3) 0.9210(8) 


z Bk 




1.426 GeV 




Vfli 1 ) 




0.8758(25) 0.8554(17) 0.8187(13) 




it 4) 




1.1865(38) 1.1496(32) 1.1241(24) 


z Bk 


44) 

4^) 
it 4) 


3GeV 


0.9765(1) 0.9549(1) 
0.9396(2) 0.9153(1) 
0.8795(4) 0.8537(2) 
1.0432(4) 1.0238(2) 



TABLE XIX. Bk renormalization factors in the four intermediate RI/SMOM schemes at the scales /I. Here 
the quoted error contains only the statistical contributions from the amputated vertices, not the fluctuations 
from the uncertainties on the lattice spacings. Note that we did not calculate the 32ID renormalization 
factors at 3 GeV as this point lies beyond the Rome-Southampton window on this lattice. 

tions for the renormalization factors that include the fluctuations on the lattice spacings, following 
the procedure in section IVI A 41 As before, separate distributions were obtained for each of the 
three chiral ansatze, with the central values shifted appropriately, allowing us to later separate the 
chiral and finite-volume systematic errors. The formation of the superjackknife distributions re- 
quires the derivatives of Zb k with respect to the lattice spacings, which we again determined by 
measuring the differences in the central values as the lattice spacings are varied by their total error. 
We use the full superjackknife distributions to renormalize Bk in the following sections. 
We determined the step-scaling factors by taking the continuum limit of the ratio of Zg K at 3 GeV 
and 1.4 GeV in each of the four schemes. The results are given in table IXXl 
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Projector P 


Chiral Ansatz 




ChPTFV ChPT Analytic 


(4,4) 


0.9140(34) 0.9145(33) 0.9150(34) 




0.9589(21) 0.9591(21) 0.9593(21) 


(4,r) 


1.0127(74) 1.0127(75) 1.0128(75) 


W4) 


0.8641(80) 0.8647(80) 0.8654(81) 



TABLE XX. Non-perturbative step-scaling factors for each intermediate scheme SMOM(P), used a poste- 
riori to run Zb k from 1.426 to 3 GeV. A different value is obtained for each determination of the lattice 
spacings. 

B. Chiral/continuum Fits 

The determination of on the 32ID ensemble set was discussed in section|Tn]and the values listed 
in tables [VTlTI and ITX1 These data and those on the Iwasaki ensemble sets were renormalized into 



Parameter 


ChPT ChPTFV 


Parameter 


Analytic 


£ 2 /dof 


0.71(45) 0.56(40) 




0.49(33) 


B 

f 


4.144(89) 4.110(93) 
0.1221(29) 0.1259(28) 






B K 

C B K ,a 

„ID 

C B K ,a 
c B K ,m x 

c B K ,m y 
c B K ,m h 


0.580(10) 0.584(10) 
0.073(44) 0.072(44) 
0.099(23) 0.095(23) 
0.00458(72) 0.00398(76) 
-0.0079(16) -0.0079(17) 
1.440(39) 1.450(40) 
-0.08(13) -0.06(13) 


/~<Bk 
/^Bfc,i 

r B K ,ID 

cf K 

r->Bfc 
c 2 

/~>Bk 
c 3 

/^Bk 
<- 4 


0.597(11) 
0.059(46) 
0.086(23) 
0.33(24) 
-0.07(54) 
1.450(40) 
-0.04(13) 



TABLE XXI. The ^ 2 /dof and parameters for each of our chiral fit ansatze for Bk , with the fits performed to 
data renormalised in the SMOM(^, d) scheme with a cut on data with corresponding pion masses m n > 350 
MeV. The parameters are given in physical units and with the heavy quark mass expansion point adjusted to 
the physical strange quark mass. For the ChPT and ChPTFV ansatze the chiral scale A^ has been adjusted 
to 1 GeV. 
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Ansatz 



Scheme 
SMOM(d,d) SMOM^,/^ 



Analytic 0.5978(87) 0.5506(77) 
0.5871(84) 0.5410(75) 
' 0.5904(85) 0.5436(75) 

TABLE XXII. Predictions for B K in the continuum limit in the SMOM(^, <f) and SMOM^ , ) schemes at 
ji = 1 .426 GeV for each global fit ansatz. These results were obtained using simultaneous/chiral continuum 
fits to renormalised data with a pion mass cut of 350 MeV. 

the RI/SMOM intermediate schemes at /i = 1.426 GeV using the results of the previous section. 
Anticipating the discussion in the following section, we present only the results of fitting to data 
renormalized in the SMOrv^y", y") and SMOM(^,^) intermediate schemes. 
As before, we obtain our chiral/continuum fit forms by performing an expansion in the quark 
masses and a 2 to NLO, with the light-quark mass expanded about both the chiral limit - using 
chiral perturbation theory - and about a fixed mass via a Taylor expansion. For example, for the 
analytic ansatz we obtain the following: 

B\ y = CS K (l + cf*' A(1) [a 1 ] 2 ) + cf **J + Cf m\ + Cf (rn) - m m ) + Cf (rn\ - m m ) , (88) 

and for the ChPT ansatz: 

R l _ pO J 1 i AW r„li2 i cbkmXi i CB K ,m x Xl Xl i ( Xx \ 1 

- a »i 1+e *>] +-^ + -^ 32^'° 8 (a|)} (89) 



+CB K ,m y (my - m/joj + c BK ,m h (m\ - mm 



1 



where % q = 2Bm q and the chiral scale A^ is set to 1 GeV. These fit forms apply specifically to 
the primary lattice 1; the forms for any other ensemble set e can be obtained by inserting factors 
of Zf and Zf and selecting the a 2 coefficient appropriate to the lattice action. The finite-volume 
correction terms for the ChPT fit form can be found by applying the rules given in Appendix C of 



ref. ni5H. 

Following the 2010 analysis strategy, we fixed the leading order LECs B and / in the ChPT fits to 
those obtained in section |V] reducing the number of free parameters. We also fix the the scaling 
factors Zj, Zf, and R a , as well as the physical quark masses and the overall scale to those obtained 
using the corresponding ansatz in section |V] 
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We once again performed cuts to the data set used in the ChPT and ChPTFV fits, reducing the 
largest pion mass to 350 MeV. In the main analysis we performed our analytic fits with a lower 
pion mass cut of 260 MeV in order to obtain a better fit to the data. When using this cut for the 
analyic fits to Bk, we found that we lost almost all statistical precision on our continuum prediction 
because the statistical errors on the 32ID ensembles become very large in the light-mass regime 
(cf. figure [20b. hence the effective number of points contributing to the fit after the cut is smaller 
than in the case of mk or fx. Raising the cut to 350 MeV produced much more reliable results, 
hence we adopt this higher cut for the analytic fits in this section. This is justified by the fact 
that we observed no statistically significant deviations of the fit functions from the data over this 
expanded range, hence we have no reason to believe that this will lead to an incorrect estimate for 
the chiral systematic error. This was not the case for the fits to m K , where we observed significant 
deviations. 

The analytic fits were again performed to data corrected to the infinite- volume using the ChPTFV 
fit form. 

The parameters and uncorrelated ^ 2 /dof obtained by fitting to data renormalized in the SMOM(^, d) 
are listed in table IXXJ and we give histograms showing the deviation of the data from the fits in 
figure[T9] We list the continuum predictions in both the SMOM(^, d) and SMOM^, y<") schemes 
in table IXXIIl 

In figure [20] we overlay the data with the fit curves on the 32ID ensembles, and in figure [2T1 we 
show the chiral extrapolation overlaying data corrected to the continuum limit via the ChPTFV 
and analytic parameterizations. The separation of the points at the physical up/down quark mass is 
used as a measure of the error on the chiral extrapolation. In these figures we see that the statistical 
errors increase substantially as we approach the chiral limit. The central values also appear to trend 
upwards, although this apparent curvature is in the opposite direction to that suggested by chiral 
perturbation theory and is therefore likely to be simply due to the low resolution on these data 
points. 

C. Final Results for Bk 

Applying the step-scaling factors given in table lXXl to the continuum predictions in table IXXlTl we 
obtained Bk in the SMOM(^) and SMOM(y M , y^ 1 ) schemes at a 3 GeV renormalization scale. 
Once again we see some cancellation between the statistical fluctuations on the step-scaling factor 
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Scheme 


Ansatz 


SMOM(^) 


SMOM(7^,7^) 


Analytic 


0.5213(72) 


0.5397(76) 


ChPT 


0.5188(72) 


0.5369(76) 


ChPTFV 


0.5282(73) 


0.5470(78) 



TABLE XXIII. Predictions for B K in the continuum limit in the SMOM(^) and SMOM^, f 1 ) schemes 
at fi = 3 GeV for each global fit ansatz. These results were obtained by applying the continuum step-scaling 
factors to the values in table IXXlTl 

and the 1 .4 GeV quantity. 

Finally, we apply the ms conversion factors given in section |VII A 21 to convert our results into the 
ms scheme for the convenience of the reader. Before quoting our final results, we first discuss the 
various contributions to the systematic error. 

1. Systematic Errors 

For our central values and statistical errors of our final ms prediction, we follow the 2010 analysis 
in taking the results obtained using the SMOM(^, d) intermediate scheme, which is best described 
by one-loop perturbation theory. Following section |V] we estimate the finite- volume and chiral 
extrapolation systematics on this quantity from the differences between the ChPTFV result (which 
we take as our central value) and the ChPT and analytic results respectively, taking for our estimate 
the larger of the superjackknife error on the difference or the difference in central values. As we 
propagated the differences between the lattice spacings through our analysis in section I VII A 4[ 
the aforementioned systematics on the renormalization factors are automatically included in the 
differences above. 

The remaining systematic errors are associated with the perturbative conversion into the ms 
scheme. The largest of these is the perturbative truncation error. To determine this we again 
follow the 2010 analysis strategy of taking the difference between the values of Bx in the ms- 
scheme at 3 GeV obtained using the SMOM(^,^) and SMOM(y M , y") intermediate schemes, 
the latter of which is also well-described by perturbation theory. As discussed in section |VI A 51 
and above, there are non-perturbative effects associated with the spontaneous chiral symmetry 
breaking and the presence of additional energy-scales (Aqcd> m s , etc.), that contribute to the per- 
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FIG. 19. Histograms of the deviation of the fit from the data for Bk on each of the three ensemble sets (321 
top, 241 middle and 32ID bottom) with the analytic (left) and ChPTFV (right) ansatze. 
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FIG. 20. The analytic (left) and ChPTFV (right) fit curves overlaying the partially-quenched data on 
the 32ID ensembles at the simulated strange quark mass. The fits were performed to the data set with 
corresponding pion masses m % < 350 MeV, with the data renormalized in the SMOM(^,^) intermediate 
scheme. 

turbative systematic. In ref. (2D we found that in the non-exceptional schemes these effects are tiny 
compared to the truncation systematic, therefore we do not include these effects in our systematic 
error budget. 



2. Final Results 

Using the ChPTFV result in the SMOM(^, d) for the central value and statistical error, and obtain- 
ing the chiral and finite-volume systematic errors as above, we find: 

5k(SMOM(^),3 GeV) =0.540(8) (7) (3). (90) 

where the errors are associated with the statistical, chiral, and finite-volume respectively. Convert- 
ing this to the Ms-scheme at 3 GeV using one-loop perturbation theory we obtain 

B k (ms,3 GeV) = 0.535(8) (7) (3) (1 1) , (91) 

where the first three errors are as before, and the final error is that associated with the truncation 
of the perturbative series. Converting to the Renormalisation-Group Invariant (RGI) scheme, we 
find 

£* = 0.758(11)(10)(4)(16). (92) 
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FIG. 21. The chiral extrapolation of Bk in the SMOM(^,^) scheme in the continuum limit. The circular 
and diamond-shaped data points in darker shades show the data corrected to the continuum limit using the 
ChPTFV fit form, and those in lighter shades via the analytic form. The circular points indicate those data 
included in the fits, and the diamond points those that were not. The upper and lower curves show the 
analytic and ChPTFV chiral fit forms and the corresponding square data points the extrapolated values at 
the physical up/down quark mass. All data and curves are shown at the physical strange quark mass. 

In the 2010 analysis we obtained: 



#k(ms,3 GeV) =0.529(5) (15) (2) (11] 



(93) 



This is highly consistent with the result of the present analysis. In our new result we see a large 
improvement in the chiral extrapolation systematic, which results from lowering the pion mass cut 
to 350 MeV from the 420 MeV used in the previous analysis. 



For comparison, the FLAG working group give Bk = 0.738(20) Il25h for Bk in the RGI scheme 
with 2+1 quark flavors, which was determined by combining our 2010 analysis result yO with the 
value calculated by Aubin et al (45J], which used domain wall valence quarks on the 2+1 flavor 
staggered fermion lattices produced by the MILC collaboration. The result of Bk = 0.758(22) 
obtained in the current analysis is consistent with this value. Other calculations performed since 
the publication of the FLAG 2010 paper include refs. D46D, D47D and 048Q. 
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Ansatz 



X 2 /dof 
Uncut 



X 2 /dof 
Cut 



Analytic 1.45(66) 0.141(71) 
1.47(67) 0.41(40) 
J 1.47(67) 0.42(40) 

TABLE XXIV. Fit ansatze and the associated uncorrelated # 2 /dof obtained by fitting to ro and r\ over the 
full data set (second column) and to the cut data set (third column). The upper bounds on the pion mass in 
the cut data sets are m K = 350 MeV for the ChPT and ChPTFV fits and m„ < 260 MeV for the analytic fit. 

VIII. CHIRAL/CONTINUUM FITS AND PHYSICAL RESULTS FOR THE SOMMER SCALES 

In this section we present the results of applying our global fit technique to the Sommer scales, ro 
and r\ . In ref. yj] we determined continuum values for these parameters using global fits to our 
Iwasaki ensemble sets. In this paper we extend these fits to include the 32ID ensemble set and 
observe the effect of lowering the pion mass cut. The values of ro and r\ measured on the 32ID 
ensemble sets can be found in section UTTl 

Assuming a linear dependence on the quark masses and on a 2 , we performed our chiral/continuum 
fits using the following form: 



r) = c ri)0 (l + c^Jfl [a 1 ] 2 ) + c rim m} + c nm {m\ - m m ) 



(94) 



on the primary lattice 1. As always the fit form describing another ensemble set, e, is obtained 
by inserting factors of Zf and Zf to convert the simulated quark masses on ensemble e into the 
matching scheme, and selecting the a 2 coefficient for the lattice action of the ensemble set. 
For convenience, we simultaneously fit both ro and r\, even though they do not share any common 
parameters other than the scaling factors, Z/ and Z/ ; . The lattice spacings and scaling factors were 
fixed to those obtained in the main analysis, with the fits repeated for each of the three chiral 
ansatze. For each fit we applied the same cuts as were performed to the data in section |Vl this 
corresponds to removing the data points on the 321, m/ = 0.008 and 241, m/ = 0.01 ensembles, 
and also in the analytic fit, the data point on the 321, m\ = 0.006 ensemble (cf. table IXlT). For 
later comparison we also quote the results of fitting to the full data set in this section, although as 
previously discussed these results are flawed due to the poor fit to several of the pion mass data 
points on the 32ID ensembles. In table IXXfVl we give the uncorrelated # 2 /dof of our fits and in 
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Uncut 






Cut 




Parameter 


Analytic 


ChPT 


ChPTFV 


Analytic 


ChPT 


ChPTFV 


c ro , (GeV- 1 ) 


2.479(34) 


2.445(36) 


2.438(38) 


2.462(49) 


2.453(51) 


2.441(52) 


c' „ (GeV 2 ) 


-0.065(53) 


-0.013(46) 


-0.008(47) 


-0.008(85) 


-0.018(64) 


-0.010(65) 


4 0)fl (GeV 2 ) 


-0.055(24) 


-0.028(26) 


-0.023(28) 


-0.032(35) 


-0.030(33) 


-0.021(34) 


c rom (GeV- 2 ) 


-1.67(87) 


-1.65(88) 


-1.64(87) 


-5.0(1.7) 


-3.6(1.4) 


-3.6(1.4) 


c, ,^ (GeV- 2 ) 


-0.83(42) 


-0.83(42) 


-0.83(42) 


-0.27(64) 


-0.56(52) 


-0.56(51) 


c n , (GeV- 1 ) 


1.697(24) 


1.675(26) 


1.671(27) 


1.662(41) 


1.650(40) 


1.642(40) 


(GeV 2 ) 


-0.099(64) 


-0.050(58) 


-0.045(58) 


0.00(11) 


0.014(91) 


0.023(92) 


(GeV 2 ) 


-0.148(25) 


-0.123(26) 


-0.118(28) 


-0.110(38) 


-0.097(38) 


-0.088(39) 


c nm (GeV- 2 ) 


-1.84(60) 


-1.82(59) 


-1.81(59) 


-2.6(2.4) 


-2.2(1.1) 


-2.2(1.1) 


c n , m , (GeV- 2 ) 


-1.02(20) 


-1.02(20) 


-1.01(20) 


-0.88(37) 


-0.73(24) 


-0.73(24) 



TABLE XXV The a 2 and mass dependences of ro and r\ obtained by fitting to the full and cut data sets. 
We repeat the fits for each choice of chiral ansatz used for the determination of the scaling parameters. The 
upper bounds on the pion mass in the cut data sets are m n = 350 MeV for the ChPT and ChPTFV fits and 
m % < 260 MeV for the analytic fit. The parameters are given in physical units and with the heavy quark 
mass expansion point adjusted to the physical strange quark mass 





Uncut 


Cut 




Analytic ChPT ChPTFV 


Analytic ChPT ChPTFV 


continuum 
'0 


2.475(33) 2.441(35) 2.435(37) 


2.451(48) 2.445(49) 2.433(50) 


^.continuum 
'1 


1.693(23) 1.671(24) 1.666(25) 


1.657(38) 1.645(38) 1.637(39) 


(n/ r o) continuum 


0.684(8) 0.684(8) 0.684(8) 


0.676(11) 0.673(11) 0.673(11) 



TABLE XXVI. Continuum predictions for r$ and r\ in GeV as well as their ratio, using scaling parameters 
obtained from each of the three global fit ansatze. The first set of columns contain the values obtained by 
fitting to the full data set, and the second set those obtained by fitting to the cut data set. The upper bounds 
on the pion mass in the cut data sets are m n = 350 MeV for the ChPT and ChPTFV fits and m n < 260 
MeV for the analytic fit. 
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figure [22] we show histograms of the deviations of the data from unity for the fits. We list the fit 
parameters in table IXX Vi and the continuum predictions for r\, tq and their ratio in table IXXVTl 
In the 2010 analysis we remarked on a tension between the fit and the value of r\ on the heaviest 
241 ensemble, which led to us inflating the error on the prediction for this quantity. In figure [23] 
we see the large apparent difference in the slopes of r\ with respect to m/ between the two Iwasaki 
ensemble sets that was responsible for this tension. It appears however that the slopes of r\ agree 
very well between the 321 and 32ID ensemble sets, which has led to a substantially better fit to 
r\ upon including the 32ID data. Restricting the fits to lighter data markedly improves our fits, 
reducing the # 2 /dof by at least a factor of three. The chiral behaviour of the data, as illustrated in 
the right-hand plots in figure [23] is now very linear. As a result of these observations, we decided 
that inflating the error on r\ is no longer necessary. 

We obtain continuum predictions for r\ , tq and their ratio from the cut fit results using the strategy 
detailed in section IVBl We find 

n = 1.637(39)(20)(8) GeV 1 = 0.3230(77) (39) (16) fm, 

r =2.433(50)(18)(13) GeV 1 = 0.4795(99) (35) (26) fm, (95) 

n/ro =0.6729(109)(30)(2), 

where the errors are statistical, chiral and finite-volume respectively. The values determined in 
ref. lit were n = 0.3333 (93) (2) ( 1 ) fm, r = 0.4870(89) (2) (2) fm and n /r = 0.6844(97) ( 1 ) (0) . 
By comparing the results in table IXXVll with those obtained in the 2010 analysis we find that, 
as with the Omega mass, the use of the generic scaling procedure for determining the scaling 
factors leads to considerably larger chiral and finite-volume systematic errors than the fixed tra- 
jectory approach. In the case of r\, we see a reduced systematic error in the continuum predic- 
tion due to the improved control over the chiral extrapolation. However for ro - which formerly 
did not display any tensions with the linear ansatz requiring error inflation - this is offset by 

laboration recently obtained 



the reduction in the amount of data. For comparison, the MILC co . 
r\ = 0.3106(17) fm \l6\ and in an earlier work ro = 0.462(12) fm [49], both of which appear to 
be consistent with our results. 
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[Y - Y^/ay 



1 2 
(Y - Y & )/ay 



FIG. 22. Histograms of the deviation of the fit from the data for rO and r\ over all three ensemble sets, 
fitting with the analytic (left) and ChPTFV (right) ansatze to the uncut (top) and cut (bottom) data sets. 
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FIG. 23. The chiral extrapolation of ro (top) and r\ (bottom) using the analytic and ChPTFV ansatze. 
The plots on the left show the fits to the full data set and those on the right to the cut data sets. We have 
overlayed the fit curves with the data points corrected to the continuum limit and physical strange quark 
mass using each of the aforementioned fit functions; those points shown in bold colors were corrected using 
the ChPTFV fits and those in pastel colors using the analytic fits. The circular data points are those included 
in the fits and the diamond points those that were not. The square points show the predicted value at the 
physical point. 
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IX. CONCLUSIONS 

Using the Iwasaki gauge action with the addition of the DSDR term we were able to simulate 
with domain wall fermions (DWF) at a relatively strong coupling (/3 = 1.75, cT l = 1.37(1) GeV) 
while retaining good chiral symmetry and topological tunneling; this enabled us to work with a 
large enough physical volume ([4.61 fm] 3 ) to accomodate pions as light as 143(1) MeV without 
suffering from large finite- volume effects {m n L ^3.2 for the lightest partially-quenched point 
and m K L rs 4 for the lightest unitary point) and without having to simulate with a large number of 
lattice sites; the dimensionless lattice volume is 32 3 x 64 x 32, where the final number is the length 
of the fifth dimension L s that governs the size of the chiral symmetry breaking in the domain wall 
formulation. 

The aim of this paper was to combine these data in a simultaneous chiral/continuum fit with our 
24 3 x 64 x 16 and 32 3 x 64 x 16 DWF ensembles with the Iwasaki gauge action at /3 =2.13 
(a~ l = 1.75(4) GeV) and /3 =2.25 (a -1 =2.31(4) GeV) respectively, and under the constraint of 
universality obtain continuumpredictions for various quantities. In this we broadly followed the 
strategy of our 2010 analysis tlLl2|]. 

The fits were performed assuming three forms for the mass dependence: the ChPTFV and ChPT 
forms were obtained from NLO SU(2) chiral perturbation theory with and without finite- volume 
corrections respectively, and the analytic ansatz from a linear Taylor expansion about an unphysi- 
cal mass point. 

The largest change from our 2010 analysis strategy was the use of the 'generic scaling' method to 
obtain the scaling parameters Z/ and Z/ 7 that relate the physical quark masses between our ensemble 
sets, and R a that relates the lattice scales. In this approach (which was discussed in ref. [lj] but 
not used in the final analysis) the scaling parameters are left as free parameters in our fits and the 
results are those that, along with the mass dependences and a 2 dependence, minimise the global 
X 2 . In the 2010 analysis we used the 'fixed trajectory' approach in which the ensemble sets were 
matched at an unphysical mass point prior to performing the fits. Changing to the generic scaling 
approach allows for differences between the scaling parameters Z/, Zj 7 and R a , which relate the 
physical quark masses and lattice spacings between the ensemble sets, as we go between the three 
chiral ansatze. We associated these with chiral and finite- volume systematic errors; in the fixed 
trajectory approach these differences would have been absorbed by other parameters in the fits. 
These differences gave rise to larger systematic errors on the lattice spacing predictions due to 
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their influence on the fit form for the Omega baryon mass, which we used to set the overall scale. 
In these fits we were able to determine the a 2 dependence of the 32ID ensembles even without a 
second lattice spacing using this action. This is because, within our power counting, the choice 
of action affects only the coefficient of the a 2 term. As all other parameters are shared with the 
Iwasaki ensemble sets, this only introduces one additional parameter per quantity. This parameter 
could in principle be determined by comparing a single data point to the continuum value predicted 
using the Iwasaki ensemble sets alone, however we choose to maximize the use of our data by 
including it in the global fit. 

We investigated removing data associated with the heavier pions, constraining our fits to a smaller 
range. For the ChPT and ChPTFV fits, we lowered the pion mass cut to 350 MeV, down from 
the 420 MeV used in the 2010 analysis. For the analytic fits, we found large deviations of our fits 
from the data when fitting to this range, necessitating a further reduction in the largest pion mass 
to 260 MeV. With this cut the analytic fit produced results with errors only slightly larger than the 
ChPT determinations. The necessity of lowering the cut for the analytic fits hints at the presence 
of non-linearity in our combined data set, which appears to be consistent with NLO SU(2) ChPT, 
although we cannot rule out other higher-order terms such as m 2 with our present statistics. 
We presented the results of simultaneously fitting m K , m^, fx, /k and in section |V] As in the 
2010 analysis, the pion, kaon and Omega baryon masses were used to set the up/down quark mass, 
strange quark mass and lattice scale respectively. We were then able to make predictions for the 
other physical quantities. For the pseudoscalar decay constants, we obtained f K = 127 '. 1 (3.8) MeV 
and fx = 152.4(3.4) MeV. These agree very well with the known continuum values of [|50D 
f n - = 130.4(2) and f K - = 156.1(8), which is a marked improvement from the 2010 analysis, 
in which the predictions for these quantities were considerably lower. The improvement stems 
mainly from our removal of data associated with the heavier pions. 

Combining our ChPT and ChPTFV fit results, we obtained values for the effective chiral couplings 
Is = 2.91(24) and I4 = 3.99(18), which we found to be highly consistent with our 2010 analysis 
results and with other lattice calculations. 

In section [VI] we discussed the renormalization of the physical quark masses into the ms scheme. 
We used variants of the Rome-Southampton RI/MOM scheme with symmetric kinematics as in- 
termediate non-perturbative schemes, which were applied at 1.4 GeV and the results run to 3 
GeV using continuum step-scaling factors. These were then converted into ms using perturba- 
tion theory. This analysis improved on the 2010 result in the use of twisted-boundary condi- 



86 



tions to remove 0(4)-breaking lattice artifacts in our measurements. We also increased the renor- 
malization scale from 2 GeV to 3 GeV, as this considerably reduces the systematic error arising 
from the truncation of the perturbative series. We obtained m M ^(Ms, 3 GeV) = 3.05(10) MeV and 
m, s (Ms, 3 GeV) = 83.5(2.0) MeV for the average up/down quark mass and strange-quark mass re- 
spectively. 

In section IVIII we applied our chiral/continuum fits to the neutral kaon mixing parameter. This 
analysis improved on the 2010 result through the inclusion of the Iwasaki+DSDR ensembles. We 
found a marked improvement in the chiral extrapolation systematic due to the inclusion of these 
data. For our final result we obtained B K (ms, 3 GeV) =0.535(16). 

Finally, in section IVITll we performed chiral/continuum fits to the Sommer scales tq and r\, for 
which we obtained r$ = 0.480(11) fm and r\ = 0.323(9) fm. Here the inclusion of the 32ID 
ensembles provided considerably greater stability to the fits than in the 2010 analysis, resulting in 
much reduced errors, particularly for r\, for which we were formerly forced to inflate the errors 
due to the poor % 2 ldof on the fits. 

Although the inclusion of the 32ID ensembles resulted in considerable improvements in the chiral 
extrapolation systematic error in most cases, there is still room for improvement. Our collaboration 
has recently gained access to IBM Blue Gene/Q computers, which have performances in the region 
of several hundred Teraflops per rack. Particularly when used with the improved techniques that 
we and others have developed (some of which are discussed in section fll] and Appendix lAl). these 
computers have the capability of generating domain wall fermion ensembles with physical quark 
masses and large enough L s and physical volumes to maintain small chiral symmetry breaking and 
finite- volume corrections. With such ensembles the necessity of extrapolating to the physical point 
will be removed and only the continuum extrapolation will remain. However, in the mean-time 
the results of this analysis, particularly the physical quark masses and lattice spacings, will be 
essential for any physics measurements performed on the Iwasaki and Iwasaki+DSDR ensembles. 
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Appendix A: Numerical Integration Scheme 

Integrators used in lattice simulation must be both reversible and symplectic. Consider a general 
Hamiltonian with both a kinetic (T) and potential (5) term: 

H=T(p)+S(U). (Al) 

In general this Hamiltonian cannot be integrated exactly, as the corresponding time evolution 
operator, 

exp ( xH] = exp ( x ( f + s) ) , (A2) 



involves non-commuting operators T and S. However, by making use of the Baker-Campbell- 
Hausdorff (BCH) formula one can separate T and S and integrate them at different steps. 
One of the simplest integrators that can be constructed in this way is the leapfrog integrator, 



UqpqW = exp Qrf ^ exp (tS) exp (\xf\ . (A3) 
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Using the BCH formula it can be shown that 

t/ QPQ (T) = exp (t (f+s) + ^(t 3 )) . (A4) 

The 0(x 3 ) error is accumulated over the integration such that the total error is &(x 2 ), hence the 
leapfrog integrator is a second-order integrator. Another popular second-order integrator is the 
Omelyan integrator, 

^qpqpq(t) = exp (axf^j exp ( ^t5J exp ^(1 - 2a) xT^j exp ( ^Sj exp (^axT^j , (A5) 



where a is a tunable parameter. 



Recent development on integrators has introduced the force gradient integrator (FGI) fl 1 811 as a 
fourth order integrator. The force gradient integrator is constructed by introducing the "force 
gradient term" into the integration steps. This extra force evaluation helps to eliminate the second 
order errors and makes the force gradient integrator a fourth order integrator. One choice of the 
force gradient integrator is 



U PGl (x) =exp | ^-^-rf j exp \^xS- l^ T 3 {5^5,r}} 

exp ( ^-xf ) exp (\xS- ^^x^S^f}}) exp ( 3 -^-xT 



(A6) 



1. Sexton- Weingarten Integration 

In practice the action contain contributions from both the gauge fields and the fermions, 

H = T(p)+S G (U)+S F (U). (A7) 

It is usually the case that the gauge force is larger than the fermion force by a factor of 10 or more. 
If both the gauge action and the fermion action are integrated in the same step then the step size 
x has to be chosen to accommodate the larger gauge force. This approach incurs an extra cost on 
the fermion part, which usually dominates the computing time. 
The Sexton-Weingarten integration scheme can be used to mitigate the issue. Define 

H=T' + S F (U) (A8) 
T'=T(p)+S G (U), (A9) 
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then T' and S-p(U) can be fit into one integrator. When integrating T' , its 2 parts T{p) and Sq(U) 
can be fit into another integrator. For example, when using the leapfrog QPQ integrator for both 
levels one has the following 

exp (rftj f^exp (l x ^'j ex P ( t ^f) ex P Q 1 ^') (A10) 
exp(Irr) »(«p(ltf)«p(ltS)e»p(i ; tf))", (All) 

where n can be chosen as any positive integer. In this way different time steps are assigned to 
Sq(U) and Sp(U), which can be tuned to minimize the cost. 



2. Hasenbusch Mass Splitting 

Hasenbusch mass splitting breaks a single fermion action into a few parts and offers a fine control 
on distributing fermion forces among them. 

The fermion action is derived from the following fermion determinant 

det ( M ^ m \ M ^) \ = r 9 ^ exp f_ t M(1) ; M t (1) ^ . (A12) 

VAft(l)Af(l) J J Y Y P V ! M^{m)M{m) K ' Y J V 



The Hasenbusch factorization II 1711 rewrites the above quotient action as a product of quotient 



actions by introducing intermediate masses 



de ' l, Mt(l)M(l) J = n de ' { Mt (m ,)M(m,) ) (A13) 
= n/*tf»* exp (-^(- ^.^J^/ '^l*) ■ (A14) 
where m = mq < m\ < ■ ■ ■ < m^+i = 1 • 

This method offers fine grained control on the sizes of the fermion forces since all intermediate 
masses miii = 1,2, • • • , k) can be tuned continuously. In what follows the symbol 5Q(m fl , m b ) will 
be used to represent the quotient fermion action 

S Q (m a ,m b ) = ^M(m b ) —r- — \— — -M 1 (m b )(j>, (A15) 

M~(m a )M(m a ) 

The Q in Sq means "quotient". Note that each quotient action has a different pseudofermion field 
0. This fact is not represented in the above symbol. 
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3. Final Scheme 

The quotient action discussed above accounts for 2 types of fermions. This is used to simulate the 
2 light quarks in our simulation. For simulating strange quark, the rational approximation needs 
to be used: 

det W^H V /2 = t 9 f 9 ) exp |%f (Mt(l)M(l)) 1/4 ! — - (M t (l)M(l)V / %^) , 

\M\\)M{\)) J V 9 y) K> ) ( M t(m)M(m)) 1/2 V K> K> ) ^ ' 

(A16) 

where rational approximations of function jc 1 / 4 and x~ 1//2 are used to evaluate the non-integer 
powers of matrices. In what follows we will use the symbol S R (m\,m2) to represent this rational 
action 

S R (mi,m 2 ) = ^ t fM t (m 2 )M(m 2 )) 1/4 ^ (M\m 2 )M(m 2 )) (A17) 

v J (Mt(mi)M(m!)) 1/2 v J 

where power functions such as x 1//4 and x -1 / 2 are understood to be shorthand notations of their 
corresponding rational approximations, the "R" in Sr means "rational". 
The final action used in the evolution contains the following components: 

H = T(p)+S G + Y,S Q (m i ^i,m i )+S R (m s , 1) +5 DSD r, (A18) 

i 

where niQ = mi, m^ + i = 1, m/ and m s represents the light quark mass and strange quark mass 
respectively. It is also possible to replace the quotient action 5q(m, 1) with two copies of the same 
rational action 5r(/«, 1). 

When evolving the above action, we use multiple levels of nested integrators to separate the differ- 
ent parts of the action. A general multi-level Sexton-Weingarten Integration scheme can be written 
as follows: 

H = T > =T > + Sl (A19) 
Ti=T{ +1 +S i+1 i = l,2,--,k-h (A20) 

where T' k = T(p). The above equations separate the entire action into k levels. 
The details of the evolution schemes for the 2 ensembles are listed in the following tables. The 
second column specifies which component of the action is used in 5,. The value given in the 
fourth column, n ; , denotes the number of integration steps for T[. This quantity is equivalent to n 
in section IaTTI 
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level(z') 


5, integrator type 


step size 


1 


5 Q (0.0042, 0.015) + S Q (0.015, 0.045) Omelyan QPQPQ 1 


1/8 


2 


5r(0.045,1)+5r(0.045,1) + S r (0.045,1) Omelyan QPQPQ 2 


3 


Sdsdr Omelyan QPQPQ 4 


4 


S G Omelyan QPQPQ 1 



TABLE XXVII. mi = 0.0042, m s = 0.045 ensemble evolution details, with total 4 levels of nested inte- 
grators. Also note that 2 copies of rational action Sr (0.045, 1) are used to replace a single quotient action 
Sq (0.045, 1). We use a = 0.22 for the Omelyan integrators. 



level(/) 


Si 


integrator type «/ 


step size 


1 


S Q (0.00 1, 0.01 )+5 Q (0.01, 0.04) 
+S Q (0.04,0.12) +S Q (0.12, 0.31) 
+S Q (0.3 1 , 0.62) + S Q (0.62, 1 ) + S R (0.045 , 1 ) 


FGI QPQPQ 3 


1/9 


2 


Sdsdr 


FGI QPQPQ 1 




3 


Sg 


FGI QPQPQ 1 





TABLE XXVIII. mi = 0.001, m s = 0.045 ensemble evolution details, with total 3 levels of nested integra- 
tors. 

Appendix B: Error Propagation in the Quark Mass Renormalization 

In section IVI A 41 we performed the continuum extrapolation of the ratios of quark mass renormal- 
ization factors, Z m [ and Z ra /,, which we defined in equation[73l These ratios combine the scaling pa- 
rameters Z/ and Z/j that represent the renormalization factors in the intermediate mass-independent 
'matching scheme' used during the fits and the non-perturbative renormalization factor Z m in the 
SMOM schemes calculated using the Rome-Southampton method. In this calculation, the propa- 
gation of statistical and systematic errors through the extrapolation and the subsequent application 
of the step-scaling factors is non-trivial. Firstly, we note that the 321 and 241 lattice spacings are 
very strongly correlated through R a (recall that a 241 is obtained as a 521 /R 2 a 41 ). As the errors on 
these quantities give rise to uncertainties on both the renormalization scale and on the coordinates 
used in the continuum extrapolation, naively treating them as independent between the lattices 
could potentially give rise to unrealistically large errors on the final renormalized quark masses. 
In the earlier parts of this analysis, the propagation of finite-volume and chiral extrapolation ef- 
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fects was performed by repeating the global fit with each of the three chiral/continuum fit ansatze 
(analytic, ChPT and ChPTFV) separately, taking the difference between these only at the final 
stage to estimate the corresponding systematic errors. All correlations were taken into account 
through the use of the superjackknife method to propagate the statistical errors. However, the 
determination of the non-perturbative renormalization coefficients was performed using bootstrap 
resampling for the error propagation. Therefore, in order to propagate the effects of the statisti- 
cal and systematic errors on the lattice spacings in a fashion consistent with the main analysis, we 
created 'superjackknife' distributions from the bootstrap distributions via the following procedure: 

1 . On each Iwasaki lattice, we calculated Zf n in each of the RI/SMOM schemes S at 1 .4 GeV 
and 3 GeV using bootstrap resampling to propagate the statistical errors. The results of 
these calculations were given in the previous section. We used only the central values of 
Zy and the lattice spacings during this procedure such that the statistical error contains only 
the fluctuations from the measurements of the amputated vertex functions. For the lattice 
spacings we used the central values from the ChPTFV determination, which we previously 
chose as our 'best' ansatz. 

2. We repeated the previous step once again, only this time we shifted the lattice spacings by 
their total error. From the change in Z s we obtained its slope with respect to a. (The slopes 
are negative for all of our schemes, as can be seen in figure [171). 

3. Using dZ^Jda we shifted Zf n to the values we would have obtained if we had repeated 
step CD using the lattice spacings obtained with the ChPT and analytic ansatze. Along with 
the original measurement we then had values of Z m with the physical scales set using the 
results of each of the three global fit ansatze. We henceforth refer to these with an additional 
superscript A denoting the chiral ansatz. 

4. For each fit ansatz we placed the corresponding bootstrap distribution on a fictitious 'su- 
perjackknife' ensemble, ensuring that the statistical fluctuations remain independent from 
others in the analysis. (Our code is able to include both bootstrap and jackknife distribu- 
tions within the same framework.) The remaining superjackknife samples were modified to 
account for the statistical fluctuations in the lattice spacings by setting each sample i to the 
following: 

(^' A ) ; = (^' A ) + ^f(^ -{«))■ 
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Here (...) denotes the central value of the distribution. 
5. For the final step we take into account the fluctuations on Zy by dividing the 'superjackknife' 

S A 

distributions for Z m ' by Zy / {Zy), where the quantity in the numerator is the superjackknife 
distribution used to normalise f n in the main analysis. 

These superjackknife distributions were used for the analysis documented in section |VI A 41 
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